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Summary  on  Studies  of  Functionally  Graded  Materials 


For  functionally  graded  materials  (FGMs),  the  change  of  microstructures 
induces  material  gradients,  and  makes  them  different  in  behavior  from  other 
materials,  homogeneous  and  conventional  composite  materials.  Our  current 
study  has  been  focused  on  the  following  two  aspects:  (1)  Micro-residual- 
stress  under  thermal  loading.  A  physical  based  computational  model  is 
developed  to  study  microstructures  in  FGMs  under  thermal  loading.  The 
influence  of  discrete  microstructure  on  residual  stress  at  the  grain  size  level 
is  examined.  (2)  Fracture  mechanics.  In  particular,  the  study  has  been  carried 
out  to  determine  the  effect  of  material  gradients  on  the  crack  tip  for  various 
loadings  and  geometries,  including  the  commonly  used  specimens,  under  the 
small-scale  yielding  condition.  A  locally  homogeneous  model  has  been 
proposed  to  predict  the  crack  propagation  direction  under  the  influence  of 
loading,  geometry  and  material  gradients.  Both  analytical  and  numerical 
methods  have  been  investigated  for  crack  mechanics  of  the  non-traditional, 
non-homogeneous  materials.  The  brief  summary  of  our  work  is  as  follows 
according  to  the  two  aspects. 


(a)  Effect  of  Microstructure  on  FGMs 

The  microstructures  of  FGMs  are  discrete  and  random  in  nature.  The 
heterogeneous  structure  causes  locally  concentrated  residual  stress  during 
thermal  loading,  as  we  have  shown  in  calculations  using  micro-mechanics 
model.  The  effect  of  discrete  microstructure  on  residual  stress  distribution  at 
the  grain  size  level  are  examined  with  respect  to  material  gradients  and  FGM 
volume  percentage  within  a  ceramic-FGM  interlayer-metal  layer  system. 
Both  thermal-elastic  and  thermal-plastic  deformations  are  considered,  and 
the  plastic  behavior  of  metal  grain  is  modeled  by  crystal  plasticity  theory. 
The  results  are  compared  with  those  obtained  using  a  continuous  model 
which  does  not  consider  the  microstructure.  In  the  averaged  sense  both 
micro-mechanics  model  and  continuous  model  give  the  same  macroscopic 
stress,  whereas  the  micro-mechanics  model  predicts  fairly  high  stress 
concentration  at  the  grain  size  level,  higher  than  700  MPa  for  300  degree 
temperature  drop  in  Ni-Al203  system.  Statistical  analysis  shows  that  the 
stress  concentration  is  insensitive  to  material  gradients  and  FGM  volume 
percentage.  This  suggests  that  the  consideration  of  microstructure  of  FGM 
for  detailed  analysis  is  needed. 


Figure  1  is  our  computational  model  in  which  the  difference  of  the 
continuous  model  and  micro-mechanical  model  is  clearly  seen.  Figure  2 
shows  the  contour  plots  of  (a)  averaged  in-plane  principal  stress  p  and  (b) 
accumulated  sum  of  slips  developed  in  the  linear  gradient,  40%  FGM.  The 
temperature  drop  was  from  700  to  400°C.  Similar  to  the  elastic  case,  the 
stress  distribution  is  inhomogeneous,  with  many  metal  grains  experiencing 
high  tensile  stresses  and  many  ceramic  grains  experiencing  compressive 
stresses.  Stress  concentration  in  the  ceramic  grains  is  significantly  reduced  in 
the  region  where  metal  content  is  greater  than  70vol%  due  to  plastic 
relaxation.  With  only  300  degree  temperature  drop,  Figure  2(b)  shows  that: 
(a)  there  is  plastic  strain  accumulation  in  many  of  metal  grains,  and  (b) 
certain  sites  have  relatively  high  strain  accumulations,  about  1.5%.  The  high 
strain  accumulation  seems  to  appear  in  the  regions  where  metal  content  is 
between  50  to  75%.  Figure  3  shows  the  distribution  profiles  of  averaged  in¬ 
plane  stress  in  FGM  layer  for  both  elastic  and  plastic  cases.  The  plastic 
relaxation  effect  is  very  clear  here  since  stresses  are  in  general  shifted  to 
lower  magnitudes.  The  distribution  for  high  tensile  stresses  with /?>700Mpa, 
however,  has  reduced  only  slightly  with  plastic  relaxation.  Similar  to  the 
thermal-elastic  case,  the  stress  distribution  profile  for  high  tensile  stress 
regions  is  insensitive  to  material  gradient  and  FGM  volume  percentage.  On 
the  other  hand,  the  distribution  profile  for  high  compressive  stresses  with 
p<-2500Mpa  (mostly  in  ceramic  grains)  drops  significantly.  This  suggests 
that,  when  ceramic  grains  are  subject  to  tensile  stresses  if  temperature 
increases,  the  plastic  relaxation  effects  may  reduce  their  tensile  stress 
concentrations.  We  average  stresses  over  each  column  of  elements  for  the 
plastic  solution  to  obtain  the  averaged  in-plane  stresses.  Compared  with  the 
elastic  averaged  in-plane  stresses,  it  is  found  that  the  metal  rich  section  and 
part  of  the  pure  metal  region  are  under  general  macroscopic  yielding,  which 
sets  the  maximum  magnitude  of  the  macroscopic  stresses  for  the  plastic 
case. 

We  also  employ  APSP  (averaged  peak  stress  of  p)  to  treat  the  data 
obtained.  The  term  6%  APSP  is  the  stress  p  averaged  over  the  6vol% 
microstructure  of  the  FGM  layer  which  has  the  highest  tensile  stresses  (p); 
similar,  one  can  obtain  3vol%  APSP  etc.  Figure  4  shows  the  6vol%  APSP 
for  different  material  gradients  and  different  FGM  volume  percentages,  and 
for  both  elastic  (marked  with  EL)  and  plastic  solution  (marked  with  PL). 
From  it,  the  distribution  profile  for  high  stresses  is  found  to  insensitive  to 
material  gradients  and  FGM  volume  percentages,  and  the  plastic  relaxation 


effect  is  relatively  small  for  high  stresses.  Similar  conclusion  can  be 
obtained  for  3vol%  APSP. 


(b)  Cracks  in  FGMs 

In  the  crack  problems  we  have  studied,  a  closed  form  solution  for  a  semi¬ 
infinite  crack  in  a  FGM  strip  is  obtained.  From  the  fundamental  solution,  we 
have  found  these  useful  solutions:  four-point-bending  specimen,  orthotropic 
cracked  FGM  plate.  These  results  show  that  the  material  gradients  do  have  a 
strong  effect  on  the  stress  intensity  factors  and  the  mode  mixity  which 
measures  the  proportion  of  mode  I  to  mode  II  at  the  crack  tip.  The 
magnitude  of  the  stress  intensity  and  the  mode  mixity  are  considered  as  the 
most  important  factors  to  determine  crack  propagation  in  FGMs.  The 
documented  complete  solution  can  be  used  in  fracture  testing  of  FGMs.  For 
the  small-scale  yielding  crack,  crack  deflection  initiation  in  FGMs  is  studied 
by  the  locally  homogeneous  model.  In  the  model,  the  effect  of 
microstructure  at  the  tip  region  is  neglected  and  as  in  all  of  our  work  for 
crack  problems  up  to  now  the  FGMs  are  considered  as  perfectly  non- 
homogeneous  materials  with  gradients  of  material  properties  at  macroscopic 
level.  The  locally  homogeneous  model  neglects  the  second  order  effects  and 
is  a  first  order  approximation.  Using  it,  we  have  examined  the  crack 
propagation  direction  in  several  cracked  specimens.  In  our  numerical  study, 
a  simplified  method  is  found  for  calculating  crack  tip  field  of  FGMs  in  finite 
element  analysis.  We  show  that  the  standard  domain  integral  is  sufficiently 
accurate  when  applied  for  FGMs  at  the  small  domain  near  crack  tip,  and  the 
non-homogeneous  term  is  very  small  compared  to  the  standard  domain 
integral.  We  have  given  the  error  estimation  in  terms  of  domain  size, 
material  properties  and  their  gradients.  The  numerical  results  for  both  two- 
dimensional  and  three-dimensional  problems  show  that  the  method  is 
accurate  and  efficient.  The  advantage  of  the  method  is  that,  it  does  not 
require  the  input  of  material  gradients  and  the  existing  finite  element  codes 
can  be  used  for  FGMs  without  much  additional  work. 

Figure  5  and  6  show  some  of  the  cracked  geometries  we  have  solved  using 
both  analytical  and  numerical  methods.  Figure  7  shows  the  solution  of  mode 
I  stress  intensity  factor  vs.  the  position  of  the  crack  tip  in  the  FGM  for  linear 
material  variation  in  the  three-layered  three-point-bending  specimen  shown 
in  Figure  6,  where  h/H  =  0.1.  The  geometry  is  the  case  that  the  interlayer  of 


FGM  is  considerably  thin  compared  to  the  two  bulk  materials.  The  solutions 
of  this  kind  for  various  h/H  and  material  variation  form  a  complete  solution 
for  the  three  point  bending  specimen.  Usually  tough  materials  such  as  metals 
have  lower  modulus  than  brittle  materials  such  as  ceramics.  From  this  figure, 
when  the  crack  travels  from  tough  side  (the  side  with  smaller  modulus)  to 
brittle  side  (the  side  with  larger  modulus)  the  crack  tip  stresses  increases. 
When  the  toughness  of  the  two  bulk  materials  are  different,  it  is  expected  to 
vary  along  the  thickness  of  the  FGM  and  can  be  written  as  r((a-H)IH)  in 
the  FGM.  Then,  for  stable  growth  in  the  FGM  we  have:  (a)  the  energy 
release  rate  is  equal  to  the  toughness,  and  (b)  the  rate  of  change  of  energy 
release  rate  with  respect  to  the  crack  length  is  less  than  that  of  the  toughness. 
For  unstable  growth,  in  (b)  the  “less”  is  replaced  by  “greater”.  Let’s  look  a 
special  case  where  the  toughness  is  constant  across  the  thickness  of  the 
FGM.  From  the  figure,  we  see  in  this  special  case  that  when  material  #2  is 
much  softer  than  material  #1,  E1iEx  «1,  the  crack  growth  is  likely  to  be 
stable.  This  is  especially  true  when  the  crack  tip  is  close  to  material  #2. 
When  material  #2  is  stiffer  than  material  #1,  the  crack  growth  is  likely  to  be 
unstable.  In  general,  if  the  toughness  varies  with  position  and  the  crack  is 
close  to  material  #2  with  E2  /£,  « 1,  it  is  likely  to  be  stable  growth  since  the 
drop  of  the  slop  in  the  figure  is  very  strong  in  the  region  and  it  is  likely  to 
overcome  the  drop  of  the  toughness.  Figure  8  is  the  case  that  hi H  -  0.5  and 
other  parameters  are  the  same  as  Figure  7.  From  the  two  figures  we  see  that 
the  trend  of  these  curves  has  a  dramatic  change  as  the  percentage  of  the 
FGM  changes,  and  this  is  especially  true  for  those  curves  with  E2/Ex<  1 .  In 
most  part  of  Figure  8  the  stress  intensity  factor  increases  as  the  crack  length 
increases.  This  means  that  if  the  increase  of  the  toughness  at  the  crack  tip  as 
the  crack  length  increases  is  not  as  fast  as  the  stress  intensity  factor,  it  is  an 
unstable  growth  for  the  crack  tip  traveling  at  most  part  of  the  FGM  and  for 
every  E2/Ex.  In  the  discussion  of  the  crack  growth,  we  have  assumed  that 
the  crack  propagates  along  the  original  direction,  since  these  are  the  cases 
where  geometry,  loading  and  material  are  symmetric  with  respect  to  the 
crack  line.  A  first  order  approximation  model,  which  is  based  on  the  local 
homogeneity,  has  been  used  by  us  to  examine  the  crack  propagating 
direction  for  several  cracked  FGM  geometries.  The  model  also  predicts  that 
crack  grows  along  its  original  direction  when  everything  is  symmetric. 


(c)  Further  Work  under  Consideration 

Our  further  work,  including  the  work  being  carried  on,  concentrates  on  the 
following  unsolved  problems  relating  to  the  understanding  of  FGMs’ 
behavior  as  well  as  the  implication  of  these  Understandings  on  the  design 
aspects  and  material  selection  for  the  Hi-Temp  purpose: 

•  Continue  our  study  in  micro-mechanical  modeling  of  residual  stress 
under  thermal  loading.  Carry  out  more  detailed  analysis  on  the  micro¬ 
stress  concentration  by  investigating  more  FGM  systems  and  more 
geometries  so  that  more  insight  can  be  obtained  regarding  to  sensitivity 
of  local  stress  concentration  to  FGM  volume  percentage  and  material 
gradients.  Morever,  include  cyclic  loading  and  combined  thermal- 
mechanical  loading.  Advanced  aspect  on  crystal  plasticity  theory  to 
account  for  the  metal  particles’  detailed  deformation  will  be  developed. 

•  Plastic  and  visco-plastic  mechanism  at  the  crack  tip  region  in  FGMs 
under  thermal,  mechanical,  and  combined  thermal  and  mechanical 
loadings,  including  detailed  characterization  of  the  crack  tip  field  under 
these  loadings  and  for  various  material  gradient  effects.  This  will  involve 
both  analytical  and  large  scale  finite  element  computation  studies.  For  the 
non-linear  finite  element  analysis  for  non-homogeneous  cracked 
materials,  the  effective  numerical  method  will  be  investigated. 

•  Develop  a  phenomenological  theory  based  on  experimental  evidence,  or 
from  micromechanics  at  a  smaller  scale  which  can  account  for  the 
interaction  of  the  two  particles  phases,  to  characterize  the  toughness  of 
FGMs;  determine  the  crack  propagation  in  terms  of  various  parameters 
including  material  gradients  and  the  crack  tip  position;  and  also 
determine  the  crack  path,  the  direction  of  deflection,  using  second  order 
and  more  accurate  theory. 


Figure  1.  Schematic  drawings  of  ceramic/FGM/metal  structures,  with  (a)  continuous 
model  and  (b)  discrete  micro-mechanics  model  using  crystal  plasticity  theory 
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Figure  2.  Contour  plots  of  (a)  hydrostatic  stress  and  (b)  accumulated  sum  of  slip 
developed  in  the  linear  gradient  40vol%  FGM. 
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Figure  3.  Hydrostatic  stress  distribution  profiles  developed  in  the  FGM  layer  (40vol% 
FGM  with  linear  gradient)  for  both  elastic  and  plastic  cases. 
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Figure  4.  The  6vol%  APSP  for  different  gradients,  volume  percentage  and  for  both  elastic 

and  plastic  solutions. 
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Figure  7.  Stress  intensity  factor  vs.  the  crack  tip  position  inside  the  FGM  interlayer  for 
the  three  point  bending  structure  in  Figure  6  with  h/H=0.1. 
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Figure  8.  Stress  intensity  factor  vs.  the  crack  tip  position  inside  the  FGM  interlayer  for 
the  three  point  bending  structure  in  Figure  6  with  h/H=0.5. 
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Abstract — Small  crack  deflection  in  brittle  functionally  graded  materials  (FGMs)  is  studied.  The 
FGMs  are  modeled  as  simply  nonhomogeneous  materials,  i.e,,  the  effect  of  microstructure  is 
neglected  and  the  material  property  variation  is  considered  to  be  continuous.  Considering  local 
homogeneity  and  the  small  scale  inelasticity  of  brittle  materials,  the  toughness  is  taken  to  be 
independent  of  direction ;  therefore,  the  crack  propagates  along  the  direction  of  maximum  energy 
release  rate,  or  the  direction  which  gives  a  vanished  mode  II  stress  intensity  factor.  Kink  directions 
for  several  specimens  which  may  be  used  to  experimentally  study  fracture  behavior  of  FGMs  are 
calculated.  It  is  shown  that  material  gradients  have  a  strong  effect  on  the  kink  direction  when  the 
crack  is  at  the  central  region  of  a  FGM,  whereas  they  have  little  effect  when  the  crack  is  close  to  the 
boundaries  of  the  FGM.  ©  1997  Elsevier  Science  Ltd. 


1.  INTRODUCTION 

In  functionally  graded  materials  (FGMs),  although  the  absence  of  sharp  interfaces  does 
largely  reduce  material  property  mismatch,  cracks  occur  when  they  are  subjected  to  external 
loadings  (Yamanouchi  et  ah ,  1990;  Holt  el  ah ,  1993).  Fractures  induced  by  these  cracks 
may  in  large  part  determine  the  overall  mechanical  and  thermal-mechanical  responses  of 
FGMs.  The  need  to  understand,  quantify  and  improve  the  toughness  of  FGMs  has  brought 
interest  in  a  fracture  mechanics  methodology  for  such  materials.  In  our  recent  work  (Gu 
and  Asaro,  1997),  stress  intensity  factors  of  several  specimens  composed  of  FGMs  were 
solved ;  the  effect  of  material  gradients  on  near-tip  fields  was  determined ;  and  possible 
fracture  criterion  was  discussed.  In  this  paper,  we  address  crack  deflection  (or  kinking)  in 
brittle  FGMs,  i.e.,  for  crack  with  arbitrary  orientation,  we  study  the  direction  of  its 
extension  when  the  critical  condition  is  met.  Here,  brittle  FGMs  are  those  having  strictly 
linear  response.  An  example  is  the  Si-C  FGM  system  in  which  both  material  phases  are 
brittle.  For  those  FGMs  made  of  metal  and  ceramic  phases,  the  present  model  gives  an 
approximate  solution  if  the  crack  is  on  the  brittle-behaved  ceramic-rich  side.  If  the  crack  is 
at  the  metal-rich  side,  its  propagation  is  primarily  via  plastic  mechanisms.  Our  study  aimed 
at  non-linear  crack  tip  behavior  is  ongoing  (Gu  and  Asaro,  1996),  and  will  be  discussed 
elsewhere. 

The  crack  deflection  model  is  developed  in  the  same  spirit  as  that  for  homogeneous 
materials  (Cotterell  and  Rice,  1980)  and  for  bimaterials  with  interface  cracks  (He  and 
Hutchinson,  1989).  The  crack  tip  stress  and  displacement  fields  of  FGMs,  as  briefly  dis¬ 
cussed  in  Section  2,  take  the  same  forms  as  those  for  homogeneous  materials.  Based  on  this 
fact,  the  asymptotic  problem,  which  has  a  homogeneous  body,  is  employed  to  study  the 
crack  tip  behavior.  Considering  the  local  homogeneity  and  small-scale  inelasticity  of  brittle 
materials  around  the  crack  tip,  the  toughness  is  taken  to  be  independent  of  direction  at  a 
fixed  point.  It  follows  that  the  crack  propagates  along  the  direction  of  maximum  energy 
release  rate  or  the  direction  in  which  the  mode  II  stress  intensity  factor  vanishes.  The  kink 
direction  is  a  function  of  the  external  loading,  the  geometry,  and  elastic  property’s  gradients 
of  a  given  specimen.  After  a  short  discussion  on  crack  deflection  model  for  homogeneous 
materials  and  for  bimaterials  in  Section  3.1  and  Section  3.2,  we  present  that  for  FGMs  in 
Section  3.3,  in  which  the  asymptotic  problem  based  on  the  K-field  and  the  directional 
independence  of  the  toughness  are  discussed  in  detail.  Kink  angles  for  four-point  bending 
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specimen,  double-cantilever  beam  and  center  cracked  plate  are  calculated  in  Section  4,  and 
the  following  qualitative  results  are  obtained  from  the  solution.  For  the  four-point  bending 
specimen,  the  crack  intends  to  grow  to  the  more  compliant  side.  For  the  double-cantilever 
beam  and  center  cracked  plate,  when  the  crack  is  at  the  middle  of  the  specimen  or  at  the 
compliant  side,  it  intends  to  grow  to  the  more  compliant  side,  whereas  when  the  crack  is  at 
the  stiff  side  it  intends  to  grow  to  the  stiffer  side.  The  material  gradients  do  have  a  strong 
effect  on  the  kink  angle  when  the  crack  is  at  the  middle  of  the  FGM ;  but  the  effect  is  small 
when  the  crack  is  close  to  the  boundaries  of  the  FGM.  We  also  investigate  crack  propagation 
in  a  compositionally  graded  interface.  It  is  found  that  of  the  two  Dundurs’  parameters,  a 
and  fi  (see  eqn  12),  the  effect  of  the  former  on  the  kink  angle  is  stronger  than  the  latter. 

2.  CRACK  TIP  FIELDS  OF  FGMs 

For  the  purpose  of  studying  crack  kinking,  the  major  results  of  the  crack  tip  fields  in 
FGMs  are  highlighted;  detailed  discussion  can  be  found  in  Gu  and  Asaro  (1997).  The 
functionally  graded  material  shown  in  Fig.  1  is  considered  as  a  nonhomogeneous  material 
whose  material  properties  vary  continuously.  Stresses  near  the  crack  tip  have  a  square-root 
singularity,  and  singular  terms  of  the  stresses  are  of  the  form : 

= - -§=<(*)+ 4=  m+ -4=4"  0) 

sjlnr  y/2nr  yJ2nr 

where  i,  j  -  1,2, 3 ;  r  and  6  are  the  polar  coordinates  shown  in  Fig.  1.  The  dimensionless 
angular  functions  ffjj(O),  6 JJ(6)  and  are  the  same  as  those  for  homogeneous  materials. 

The  result  is  independent  of  the  form  for  material  properties  and  the  orientation  of  the 
crack.  The  stress  intensity  factors  Kh  K„  and  K„,  are  functions  of  the  material  gradients, 
external  loading,  and  geometry.  Material  gradients  do  not  affect  the  order  of  the  singularity 
and  the  angular  functions,  but  do  affect  the  stress  intensity  factors.  As  a  result  the  near-tip 
stresses  have  the  same  form  as  that  for  a  homogeneous  material.  It  can  also  be  shown  that 
the  near-tip  displacements  take  the  same  form  as  that  for  homogeneous  materials,  and  this 
is  independent  of  material  gradients  and  the  orientation  of  the  crack.  For  plane  stress  and 
plane  strain  problems,  they  are  of  the  form : 


K, 

2£'(0) 


K„ 

2£'(0) 


(2) 


where  £'(0)  is  the  Young’s  modulus  at  the  crack  tip,  and  the  angular  functions,  14(6)  and 
wf(0),  are  the  same  as  those  for  homogeneous  materials. 

Having  the  near  tip  stress  and  displacement  fields,  the  energy  release  rate  of  the  crack 
tip  is  obtained  as 


K)  ,  Kj, 


+  • 


K}„ 


£'(  0)  £'(0)  2/i(0)' 


(3) 


Here,  /t(0)  is  the  shear  modulus  at  the  crack  tip.  The  above  equation  is  again  independent 


Fig.  ].  A  crack  in  a  FGM  which  has  continuous  variation  of  material  properties. 
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of  the  form  for  the  material  properties  and  the  orientation  of  the  crack,  and  has  the  same 
form  as  that  for  homogeneous  materials.  The  path-independence  of  the  J  integral  (Rice, 
1968)  holds  if  the  crack  is  perpendicular  to  the  direction  along  which  material  properties 
change ;  this  is  implied  in  Rice’s  original  proof  for  homogeneous  materials.  Using  the  near 
tip  fields  obtained  above,  it  can  be  shown  that  the  J  integral  is  equal  to  the  energy  release 
rate  for  the  crack  perpendicular  to  the  direction  along  which  material  properties  change. 

For  in-plane  problems,  the  complex  stress  intensity  K  =  K,+  iKn  for  FGMs  has  the 
following  generic  form 


K=\K\e*, 


(4) 


where 


^  =  tan  '  — 


(5) 


is  the  phase  angle  of  the  complex  stress  intensity  factor.  The  phase  angle  measures  mode 
mixity,  i.e.,  the  proportion  of  the  shear  traction  to  the  normal  traction  ahead  of  the  crack 
tip,  since 


tan-'fM  .  (6) 

\ayy/e=0,r->0 

The  complex  stress  intensity  has  the  following  dimensional  form : 

K=  TLmY ,  (7) 

where  T  is  a  representative  stress  magnitude,  L  is  a  characteristic  length  and  Y  is  a 
dimensionless  function  which  relates  to  the  geometry  of  the  problem  and  material  proper¬ 
ties.  Both  the  phase  angle  and  the  dimensional  form  are  consistent  with  those  for  homo¬ 
geneous  materials. 


3.  MODELING  OF  CRACK  DEFLECTION 

A  number  of  studies  have  been  performed  for  crack  deflection  in  homogeneous 
materials  and  bimaterials  with  interface  cracks.  Before  discussing  the  model  for  FGMs,  we 
give  a  brief  review  of  some  results  for  homogeneous  materials  and  bimaterials  in  Section 
3.1  and  3.2,  respectively.  Since  the  interest  is  crack  extension  initiation,  we  concentrate  on 
small  kinks,  where  the  kink  length  is  much  less  than  the  length  of  the  pre-existing  cracks 
(main  cracks). 

3.1.  Homogeneous  materials 

Stress  intensity  factors  at  the  kinked  crack  tip  uniquely  characterize  its  near-tip  fields ; 
therefore,  they  are  the  parameters  to  determine  the  deflection  direction  at  the  load  when  the 
crack  starts  to  propagate.  The  stress  intensity  factors  for  the  kinked  crack  in  a  homogeneous 
material  have  been  solved  for  elsewhere,  including  Palaniswamy  and  Knauss  (1978)  and 
Lo  (1978).  The  latter  detailed  a  method  for  solving  crack  kinking  problems  using  integral 
equations  which  are  formulated  by  the  continuous  distribution  of  dislocations :  a  method 
valid  for  both  finite  and  infinitesimal  (small)  kinks.  The  idea  is  essentially  to  remove  the 
tractions  which  are  caused  by  the  stress  field  of  the  main  crack  in  the  kinked  crack  before 
kinking.  For  the  infinitesimal  kink,  since  the  kink  length  is  much  smaller  than  the  size  of 
the  K-dominance  zone  of  the  main  crack,  the  loading  is  the  two  stress  intensity  factors,  K} 
and  KJh  of  the  main  crack,  and  the  stress  intensity  factors  of  the  kinked  crack,  Kfand  Kfi , 
are  obtained  in  terms  of  K}  and  KIh  This  is  the  asymptotic  problem:  a  semi-infinite  crack 
in  the  homogeneous  material  which  is  loaded  by  the  K-field  of  the  main  crack  characterized 
by  Kj  and  Ku\  the  propagation  of  the  main  crack  is  controlled  by  its  K-field.  At  large 
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distances  from  the  semi-infinite  crack  tip,  the  stress  field  approaches  the  K-field  of  the  main 
crack.  Near  the  semi-infinite  crack  tip,  the  stress  field  is  perturbed  from  the  K-field  of  the 
main  crack  because  of  the  kinking.  By  dimensional  analysis  and  linearity,  the  stress  intensity 
factors  of  the  kinked  crack  are  expressed  as 


Kf=  Cu(<p)K,+Cl2((t>)K,„ 

KZ^C^Kj+CniW,,,  (8) 


where  <p  is  the  angle  between  the  kink  direction  and  the  main  crack,  and  Cn,  C)2,  C2t  and 
C22  are  coefficients  which  can  be  determined  by  Lo’s  method.  For  finite  kinks,  ATfand  Kf, 
must  be  obtained  by  solving  a  full  boundary  value  problem  considering  the  load,  the 
geometry  of  the  specimen  including  both  the  main  crack  length  and  the  kink  length. 

If  there  is  only  mode  I  loading,  the  crack  would  extend  along  the  direction  of  the  pre¬ 
existing  crack.  This  direction  is  the  direction  of  maximum  energy  release  rate.  For  mixed 
mode  problems,  the  two  often  used  criteria  are  the  maximum  energy  release  rate  criterion 
(Cotterell,  1965)  and  mode  I  type  criterion  (also  referred  to  as  local  symmetry  criterion, 
see  Cotterell  and  Rice  (1980)  and  Goldstein  and  Salganik  (1974)).  The  former  states  that 
the  crack  propagates  along  the  direction  of  maximum  energy  release  rate,  and  the  latter 
that  the  crack  grows  along  the  direction  for  which  the  mode  II  stress  intensity  factor 
vanishes.  Kink  directions  determined  by  the  two  criteria  are  consistent :  it  was  shown  that 
the  difference  is  less  than  1  degree  for  almost  all  loading  combinations  except  the  case  in 
which  the  shear  mode  is  overwhelmingly  dominant  where  the  difference  is  then  about  2 
degrees  (He  and  Hutchinson,  1989). 

For  small  <j>,  using  first  order  approximation,  Cotterell  and  Rice  (1980)  were  able  to 
analytically  evaluate  those  coefficients  in  (8)  as 


1  /  6  3  <b\ 

c„  =-|3C0S-+C0Syj, 

„  3  (  .  4>  .3  <p\ 

Ci  2  =  -  -lsm-+siny  1, 

_  1  /  .  4>  .  3<£\ 

C21  =-l  sm-  +  smy  1, 

1  (6  3  d>\ 

c  22  =^lcOS-+3cOSy  1. 


(9) 


They  showed  that  stress  intensity  factors  calculated  by  using  (9)  are  in  good  agreement 
with  those  exact  solutions  for  <t>  up  to  40  degrees:  the  error  is  less  than  5%.  They  also 
showed  that,  by  substituting  the  approximation  (9)  into  (8),  the  energy  release  rate  is  locally 
a  maximum  for  the  mode  I  path  (the  path  with  a  vanished  mode  II  stress  intensity  factor). 

3.2.  Bimaterials  with  interface  cracks 

For  an  interface  crack,  stresses  have  an  oscillatory  singularity,  and  both  stress  intensity 
factors  and  angular  functions  involve  Dundurs’  parameters,  i.e., 


Re(A>,£) 
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where  K  =  K,+iK„,  and 
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In  (10),  /?  is  one  of  the  two  Dundurs’  parameters.  The  Dundurs’  parameters  (Dundurs, 
1969),  a  and  /?,  are  defined  as 


Mi(x2  +  1)-^(xi  +  1) 

^,(k2  +  1)+^(k:,  +  1)’ 

~ 1)  — —  1) 
/il(K2  +  l)+//2(>M  +  1)  ’ 


where  /z,  and  /z2  are  the  shear  moduli  of  the  two  bulk  materials ;  k,  =  3 — 4v,  for  plane  strain 
and  k,  =(3-v,)/(l  +v,)  for  plane  stress  (i  =  1,  2),  with  v,  and  v2  being  the  Poisson’s  ratios 
of  the  two  bulk  materials.  The  complex  stress  intensity  factor,  K,  has  the  dimensional  form 


K=TL'l2-iEY,  (13) 

where  T  is  a  representative  stress  magnitude,  L  is  a  characteristic  length  and  y  is  a 
dimensionless  function  which  relates  to  the  geometry  of  the  problem  and  Dundurs’ 
parameters. 

For  the  interface  crack  kinking  problem,  the  stress  intensity  factors  of  the  kinked  crack 
are  not  unique  in  the  sense  that  they  are  dependent  on  the  kink  length.  For  small  kinks, 
dimensional  analysis  and  linearity  give  the  relationship  between  the  complex  stress  intensity 
factor  of  the  kinked  crack  K*  and  that  of  the  interface  crack  K  as 


K*  =  c(4>,  a,  P)Kaic  +  d(<t>,  a,  (14) 

where  a  is  the  kink  length  and  (“)  denotes  the  conjugate  of  the  complex  variable.  It  is  seen 
that  only  when  /J  =  0  (i.e.  e  —  0)  K*  is  independent  of  the  kink  length.  A  complete  solution 
in  this  case  was  obtained  by  He  and  Hutchinson  (1989)  using  the  integral  equation  method. 
The  kink  length  dependence  case,  /?  ^  0,  was  studied  by  Mukai  et  al.  (1990)  and  Geubelle 
and  Knauss  (1994).  The  results  showed  that  K*  and  the  kink  angle  are  strongly  dependent 
on  kink  length  for  sufficiently  large  /?.  It  was  concluded  by  the  latter  that  this  is  due  to 
“rotational  stress  and  deformation  fields”  (Symington,  1987)  at  the  crack  tip,  which 
extended  to  a  region  far  outside  of  the  contact  zone.  They  suggested  that  the  kink  length  a 
should  be  viewed  as  a  property  of  the  bimaterial  combination,  and  be  determined  by  fitting 
experimental  data. 

Whether  the  interface  crack  stays  on  the  interface  or  kinks  into  one  of  the  bulk 
materials  is  decided  by  the  ratio  of  the  energy  release  rate  for  the  crack  to  extend  on  the 
interface  to  that  for  the  crack  to  kink  IS s,  and  the  ratio  of  the  interfacial  toughness  T,  to 
the  toughness  of  the  bulk  materials  T,.  Kinking  is  thereby  favored  if 


(15) 


3.3.  Functionally  graded  materials 

Crack  propagation  is  the  competition  between  the  energy  release  rate  and  the  toughness 
of  the  material.  In  order  to  address  crack  kinking  in  FGMs,  we  need  to  study  both  physical 
parameters. 

The  energy  release  rate  of  a  kinked  crack  in  a  FGM  depends  on  the  geometry,  the 
loading,  and  the  material  properties  including  material  gradients.  The  exact  solution  for  it 
has  to  be  obtained  from  the  solution  of  the  full  boundary  value  problem.  For  small  kinks, 
the  solution  may  be  obtained  in  an  asymptotic  way  similar  to  that  outlined  in  Section  3.1. 
In  Section  2,  we  have  noted  that  crack  tip  fields  for  FGMs  have  the  same  forms  as  those 
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Fig.  2.  Local  homogenization  near  the  crack  tip  region  in  a  FGM :  a  homogeneous  body  is  loaded 
far  away  from  the  tip  of  a  semi-infinite  crack  by  the  K-field  of  the  crack  in  the  FGM  (the  asymptotic 

problem). 


for  homogeneous  materials,  and  only  the  stress  intensity  factors  are  functions  of  material 
gradients.  The  elastic  constants  appeared  in  the  displacement  near*tip  field  and  the  energy 
release  rate  are  those  at  the  crack  tip.  This  suggests  that,  to  study  crack  kinking  one  can 
consider  the  asymptotic  problem  shown  in  Fig.  2 :  a  homogeneous  plate,  which  has  a  semi¬ 
infinite  crack  and  the  elastic  constants  at  the  main  crack  tip,  is  subjected  to  the  loading 
characterized  by  the  stress  intensity  factors  of  the  main  crack.  The  only  difference  between 
the  asymptotic  problem  for  FGMs  and  that  for  homogeneous  materials  is  that  here  the 
elastic  constants  appearing  in  the  problem  are  those  at  the  main  crack  tip  which  change 
with  its  position,  whereas  for  homogeneous  materials  those  elastic  constants  do  not  change 
with  the  position  of  the  main  crack  tip.  Knowing  the  asymptotic  problem,  stress  intensity 
factors  of  the  kinked  crack  can  be  solved  by  the  same  technique  to  solve  those  for  homo¬ 
geneous  materials.  The  local  homogenization  results  in  that  the  relationship  (8)  for  homo¬ 
geneous  materials  holds  for  FGMs.  Specifically,  the  coefficients  C0  in  (8)  are  the  same  as 
those  for  homogeneous  materials,  and  material  gradients  affect  the  stress  intensity  factors 
of  the  kinked  crack  K*  only  through  the  stress  intensity  factors  of  the  main  crack  K.  These 
coefficients  can  be  obtained  accurately  by  using  the  integral  equation  method  to  solve  the 
asymptotic  problem  in  Fig.  2,  which  is  a  homogeneous  crack-kinking  problem ;  as  men¬ 
tioned  before,  they  are  well  approximated  by  the  expressions  in  (9)  for  small  <f>. 

The  above  approach  for  crack  kinking  is  valid  if  the  kink  length  is  sufficiently  smaller 
than  the  size  of  K-dominance  zone.  To  examine  the  K-dominance  zone,  the  asymptotic 
solution  over  the  full  field  solution  ahead  of  the  crack  tip,  a°yy ,  is  plotted  in  Fig.  3  for 
the  center  cracked  plate  subjected  to  remote  stress  shown  in  Fig.  4(d),  where  h/H  =  1.  It  is 
assumed  that  the  width  of  the  plate  is  much  larger  than  the  crack  length.  The  crack  is 
perpendicular  to  the  direction  of  material  property  variation,  and  the  variation  is  in  the 
exponential  form  which  is  the  same  as  that  in  Gu  and  Asaro  (1997)  and  which  will  also  be 
stated  in  the  next  section.  It  is  seen  that  the  size  of  the  dominance  zone  decreases  as  y, 
which  is  the  measure  of  material  gradients,  increases.  At  10%  of  the  crack  length  ahead  of 
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Fig.  3.  Comparison  of  asymptotic  solution  with  full  field  solution. 
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(a) 


(b) 


a  |  |  | 


i  i  1 

(d) 


Fig.  4.  Several  Specimens,  (a)  Three-point  bending  specimen,  (b)  Double-cantilever  beam,  (c)  Four- 
point  bending  specimen,  (d)  Center  cracked  plate.  The  material  properties  vary  along  the  vertical 

direction. 


the  crack  tip  the  error  between  the  exact  and  asymptotic  solutions  for  homogeneous 
materials  is  5%,  whereas  it  is  7%  for  yh  =  1 .2,  and  14%  for  yh  =  2.4.  Note  that,  in  the  case 
of  yh  =  1.2,  the  ratio  of  the  Young’s  modulus  at  the  upper  boundary  to  that  at  the  lower 
boundary  is  about  1 1.  It  is  fair  to  say  that,  for  material  gradients  which  are  not  too  large, 
the  radius  of  K-dominance  zone  would  be  in  the  size  comparable  to  that  of  homogeneous 
materials,  at  least,  not  reduced  much  from  that  of  homogeneous  materials.  For  interface 
cracks  in  bimaterials  and  sandwich  structures,  the  size  of  the  dominance  zone  was  inves¬ 
tigated  by  O’Dowd  et  al  (1992),  Shih  (1991)  and  Gu  (1993). 

The  often  used  techniques  to  make  FGMs  are  thermal  spray,  powder  processing  and 
chemical  vapor  deposition  (CVD).  The  microstructure  of  these  FGMs  depends  on  these 
manufacturing  processes  (see  Yamanouchi  et  al. ,  1990 ;  Holt  et  al. ,  1993).  For  a  real  FGM, 
typical  micrograph  shows  discrete  structure.  If  the  FGM  is  made  of  material  phases  A  and 
B,  the  A-rich  side  has  a  dispersive  structure  with  B  particles  in  the  A  matrix ;  at  the  B-rich 
side  A  particles  are  in  the  B  matrix ;  and  in  the  middle  region  between  the  two  sides,  it  is  a 
skeletal  structure  due  to  the  connectivity  of  both  phases.  For  such  complicated  structure, 
the  characterization  of  its  toughness  is  an  open  issue  at  this  moment.  The  toughness  is 
likely  a  function  of  the  position  of  the  crack  tip,  and  may  also  depend  on  the  direction 
along  which  the  crack  propagates  and  the  loading  phase  angle.  In  this  study,  we  neglect  the 
effect  of  microstructure ;  we  study  the  idealized  case,  simply  nonhomogeneous  materials, 
i.e.,  the  materials  are  those  with  a  continuous  change  of  material  properties.  This  means 
that  the  microstructure  of  a  FGM  is  sufficiently  fine  that  the  continuum  model  gives 
satisfactory  predictions. 

For  the  simply  nonhomogeneous  material,  since  we  study  small  kink  in  the  locally 
homogeneous  body  controlled  by  the  K-field  (the  asymptotic  problem  shown  in  Fig.  2)  and 
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small-scale  inelastic  deformation  around  the  crack  tip,  the  toughness  of  the  non- 
homogeneous  material  is  taken  to  be  independent  of  direction  at  a  fixed  point.  If  the 
cohesive  stress  p ,  the  inelastic  stress  in  the  cohesive  zone  ahead  of  the  crack  tip  to  restrain 
separation  for  creating  free  surfaces  (Dugdale.  1 960 ;  Barenblatt.  1 962)  is  not  only  a  function 
of  separation  but  also  position,  say  p(S,y),  the  toughness  of  the  FGM  at  the  position  y 
may  be  expressed  as 


r  = 


06) 


where  <50  is  limit  separation  at  y.  This  is  the  result  of  applying  the  J  integral  around  the 
cohesive  zone  of  the  homogeneous  body  in  Fig.  2  (see  Rice,  1968). 

For  some  FGMs,  their  micrograph  may  show  layered  structures,  and  these  discrete 
microlayers  have  varying  compositions  and  thus  form  the  macroscopic  material  gradients 
along  the  thickness  direction.  If  the  crack  lies  on  one  of  the  interfaces,  it  may  behave  like  a 
real  interface  crack.  The  mechanics  of  elastic  interface  fracture  was  established  (Rice,  1988 ; 
Hutchinson,  1990).  The  criterion  for  an  interface  crack  to  grow  on  the  interface  is 

^  =  ro/o,  (17) 

where  \j/  is  the  phase  angle.  The  condition  for  crack  kinking  into  one  of  the  two  adjacent 
layers  was  stated  in  Section  3.2.  Plastic  interface  fracture  was  investigated  by  Shih  and 
Asaro  (1988,  1989)  and  Shih  (1991).  In  this  study,  since  the  FGM  is  considered  as  a  simply 
nonhomogeneous  material,  the  possibility  of  having  discrete  microlayers  is  excluded. 

Considering  the  above,  the  kink  direction  for  the  FGM  is  the  direction  of  maximum 
energy  release  rate  or  that  in  which  the  mode  II  stress  intensity  factor  vanishes.  Both 
directions  relate  to  the  geometry,  loading  and  the  material  gradients.  The  two  criteria  for 
FGMs,  like  those  for  homogeneous  materials,  are  also  consistent,  since  they  are  built  on 
locally  homogeneous  materials. 

The  locally  homogenized  model  is  expected  to  work  well  for  brittle,  simply  non¬ 
homogeneous  materials,  as  discussed  above.  For  those  FGMs  in  which  plastic  mechanisms 
are  involved  in  a  sufficiently  large  region  around  the  main  crack  tip,  or  microstructural 
gradients  are  presented  at  the  tip  region,  further  investigation  is  needed  to  determine  these 
effects. 


4.  SOLUTIONS  AND  IMPLICATIONS 

An  immediate  consequence  of  the  present  model  is  that,  for  the  three  point  bending 
specimen  shown  in  Fig.  4(a),  the  crack  extends  vertically  ahead  of  the  main  crack  tip  along 
its  direction,  since  the  crack  tip  only  has  a  mode  I  stress  intensity  factor  (both  geometry 
and  loading  are  symmetric).  The  conclusion  is  independent  of  the  form  of  the  material 
property  variation,  and  is  the  same  as  that  in  the  case  of  a  homogeneous  material. 

Stress  intensity  factors  for  the  double-cantilever  beam,  Fig.  4(b),  and  the  four-point 
bending  specimen.  Fig.  4(c),  were  obtained  by  Gu  and  Asaro  (1997),  where  material 
properties  were  assumed  to  follow  the  exponential  form : 

E'(y)  =  E0  ev-v, 

v'(y)  =  VoO+wOe75’.  (18) 

In  (18),  y  and  p  are  material  constants  representing  the  material  gradients;  £0  and  v0  are 
the  values  of  these  elastic  properties  at  y  =  0.  For  plane  stress  problems,  E\y)  «  E{y)  and 
v'OO  =  v(v),  where  E(y)  and  v(v)  are  Young’s  modulus  and  Poisson’s  ratio,  respectively; 
for  plane  strain  problems,  E'(y)  =  £(v)/[l  —  v(v)2]  and  v'(y)  =  v( v)/[l  -v(y)].  The  shear 
modulus,  p(y),  relates  to  Young’s  modulus  and  Poisson’s  ratio  by 
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Ky) 


£iy) 

2[l+v'0’)]' 


(19) 


The  above  forms  provide  analytical  flexibility  and  lead  to  somewhat  simple  forms  for  the 
field  equations.  Using  (18)  and  (19),  it  is  shown  that  for  a  traction  boundary  value  problem, 
the  stress  field  depends  on  the  material  parameter  y,  which  is  related  to  the  moduli  at  the 
upper  boundary  E'u  and  at  the  lower  boundary  E\  as 


/  h ,  £“ 

zV 


In  (20),  L  =  h  +  H  is  the  thickness  of  the  FGM  and  h  a  characteristic  length.  The  parameter 
p  in  (18)  does  not  affect  the  solution.  When  the  modulus  at  the  upper  boundary  is  larger 
than  that  at  the  lower  boundary,  for  example,  the  upper  side  is  ceramic  and  the  lower  side 
is  metal,  the  parameter  yh  is  larger  than  zero ;  for  the  case  of  a  homogeneous  material,  it  is 
zero.  We  examine  the  kink  direction  in  the  double-cantilever  beam  and  the  four-point 
bending  specimen. 

For  the  double-cantilever  beam,  the  energy  release  rate  normalized  by  M/hhS,  is 
plotted  vs  the  possible  kink  angle  $  in  Fig.  5,  where  (f)  is  positive  when  the  crack  goes 
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Fig.  6.  Kink  angle  for  the  double-cantilever  beam. 


downwards  (Fig.  2).  When  h/H  =  1  (the  crack  is  at  the  middle  of  the  strip)  and  h/H  =  10 
(the  crack  is  at  the  compliant  side),  the  energy  release  rate  is  enhanced  if  (/>  >  0;  and  it 
increases  as  the  material  gradient,  measured  by  EJE'h  increases.  When  h/H  =  0.1  (the  crack 
is  at  the  stiff  side),  the  trend  is  different :  the  energy  release  rate  is  enhanced  if  (j>  <  0 ;  it 
decreases  as  the  material  gradient,  E'JE',,  increases.  For  some  crack  positions,  the  material 
gradient  does  have  a  strong  effect  on  the  energy  release  rate.  For  example,  at  h/H  =  1,  the 
maximum  energy  release  rate  for  E'JE'i  =  20  is  about  20,  whereas  it  is  12  for  the  homo¬ 
geneous  case,  E'JE]  =  1 .  The  implication  of  these  qualitative  behaviors  is  that  when  the 
crack  is  at  the  position  where  h/H  is  sufficiently  large,  the  crack  is  kinked  to  the  more 
compliant  side,  or  downwards ;  when  the  crack  is  at  the  position  where  h/H  is  small,  the 
crack  is  kinked  to  the  stiffer  side,  or  upwards.  For  the  case  of  h/H  =  1 ,  the  kink  angle,  the 
angle  at  which  maximum  energy  release  rate  occurs,  is  0°  for  homogeneous  materials.  This 
is  a  well  known  result,  and  is  due  to  the  symmetric  loading  and  symmetric  material 
properties  in  this  case.  It  is  also  seen  that,  for  the  cases  of  h/H  =  1  and  h/H  =  10,  the  kink 
angle  increases  as  the  material  gradient  increases.  Kink  angles  for  the  three  crack  positions 
are  plotted  in  Fig.  6.  It  shows  that  for  small  y,  the  kink  angles  for  the  three  crack  positions 
are  quite  different,  but  as  y  increases  they  become  closer.  When  h/H  =  0.1,  the  kink  angle 
changes  from  negative  sign  to  positive  sign,  i.e.,  upward  kink  becomes  downward  kink,  at 
yh  %  1.5  as  it  increases.  But  at  this  crack  position,  it  is  seen  from  (20)  that  E'JE',  =  40  for 
yh  =  0.34.  This  tells  us  that  the  portion  of  the  curve  where  yh  >  0.34  is  no  practical  use, 
since  E'JE',  >  40  in  that  portion  and  such  large  difference  in  elastic  moduli  may  not  appear 
in  a  FGM,  at  least,  at  the  present  time.  It  is  noted  that  when  the  crack  is  at  the  middle 
position,  the  material  gradient  does  have  a  strong  effect  on  the  kink  angle,  whereas  it  has 
little  effect  for  the  crack  close  to  the  upper  and  lower  boundaries. 

For  the  four-point  bending  specimen,  it  is  assumed  that  the  horizontal  distance  between 
the  crack  tip  and  the  loading  is  large  enough.  As  shown  in  Fig.  7,  the  maximum  energy 
release  rate  is  attained  when  <j>  is  positive,  and  increases  as  the  material  gradient,  E'JE',, 
increases  for  h/H  —  1,  0.1  and  10.  When  h/H  =  0.1  and  1,  the  increase  of  it  is  significant. 
For  example,  at  h/H  =0.1,  the  energy  release  rate  is  0.003  for  the  case  of  homogeneous 
materials,  and  0.017  for  E'JE',  =  20;  at  h/H  =  1,  it  is  8.41  and  20.02  for  the  two  cases, 
respectively.  Kink  angles  for  the  three  crack  positions  are  plotted  vs  the  material  gradient 
yh  in  Fig.  8.  The  figure  shows,  for  the  four-point  bending  specimen,  yh  does  not  have  a 
strong  effect  on  the  kink  angle.  When  the  crack  is  at  the  middle  of  the  beam,  the  difference 
of  it  between  the  case  of  a  homogeneous  material  and  that  of  E'JE',  =  20  is  less  than  2°. 
When  the  crack  is  at  h/H  =  0.1  and  10  the  difference  is  even  less.  In  experiments,  such  small 
difference  is  difficult  to  be  measured. 

For  the  center  cracked  plate  shown  in  Fig.  4(d),  the  energy  release  rate,  besides  h/H 
and  yh,  also  varies  with  the  crack  length  a/h.  For  the  case  of  a/h  =  1 ,  kink  angles  are  shown 
in  Fig.  9.  The  trend  of  these  curves  is  similar  to  that  in  Fig.  6  for  the  double-cantilever 
beam. 
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Fig.  8.  Kink  angle  for  the  four-point  bending  specimen. 


Two  layers  can  be  jointed  by  a  compositionally  graded  material  (FGM)  sandwiched 
between  them,  as  shown  in  Fig.  10.  It  reduces  the  mismatch  of  the  bimaterial  without  the 
interlayer  and  gives  better  thermal-mechanical  performance  for  the  whole  system  (Gian- 
nakopoulos  et  al .,  1995).  Also,  for  a  bimaterial,  the  interface,  the  transition  zone,  is  a  FGM 
microscopically.  For  some  bimaterials,  the  thickness  of  the  zone  is  at  the  level  that  it  may 
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Fig.  9.  JCink  angle  for  the  center  cracked  plate. 


not  be  neglected  in  analyzing  the  fracture  of  the  interface.  We  use  the  above  method  to 
study  the  crack  propagation  inside  the  compositionally  graded  layer.  The  main  crack  shown 
in  the  above  figure  is  parallel  to  the  boundaries  of  the  graded  layer.  The  thicknesses  of  the 
two  bulk  materials  are  considered  as  much  larger  than  the  thickness  of  the  interlayer.  For 
a  crack  on  a  sharp  interface,  as  discussed  in  Section  3.2,  the  energy  release  rate  of  the 
kinked  crack  is  dependent  on  the  extension  length,  and  there  is  a  contact  zone  in  the  region 
very  close  to  the  main  crack  tip  due  to  its  oscillatory  singular  field.  The  solution  to  the 
crack  deflection  problem  in  such  case  is  not  unique,  except  one  of  the  two  Dundurs’ 
parameters,  /?,  vanishes.  It  is  seen  that  by  introducing  the  transition  zone  to  the  interface 
the  oscillatory  singular  field  near  the  main  crack  tip,  which  is  physically  unacceptable,  is 
removed,  and  therefore  the  problem  of  non-uniqueness  for  crack  kinking  does  not  exist. 

The  thickness  of  the  graded  layer  and  the  position  of  the  crack  inside  it  are  specified 
by  h  and  H  shown  in  Fig.  10.  Here,  since  the  thickness  of  the  graded  layer  is  considered 
much  less  than  the  length  scale  of  the  two  bulk  materials,  the  loading  is  represented  by  the 
complex  stress  intensity  factor  K *  of  the  bimaterial  problem  without  the  graded  layer.  The 
complex  stress  intensity  factor  K  of  the  crack  inside  the  graded  layer  is  expressed  as 

K  =  qei'0hieKx.  (21) 

The  equation  is  obtained  by  dimension  analysis  and  linearity  consideration.  Both  q  and  co 
in  above  expression  are  functions  of  Dundurs’  parameters,  a  and  /?,  the  position  of  the 
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crack  inside  the  graded  layer,  hjH ,  and  the  form  of  the  material  gradients.  The  former  can 
be  evaluated  by  using  the  J  integral.  For  the  case  that  the  elastic  modulus  follows  the 
exponential  form  given  in  (18)  and  Poisson’s  ratio  is  a  constant, 

q  = - yfidL - - - .  (22) 

For  h\H  =  1,  Yang  and  Shih  (1994)  gave  a  general  expression  to  estimate  co  for  any  form 
of  material  property  variation  in  terms  of  a  known  bimaterial  solution.  They  compared  the 
results  obtained  from  the  approximation  with  those  from  finite  element  analyses,  and  it 
was  concluded  that  the  approximation  was  quite  satisfactory.  We  use  their  expression  to 
examine  crack  deflection  in  the  graded  layer  for  the  case  where  hjH  =  1  and  material 
properties  follow  the  above  mentioned  form.  In  Fig.  11,  we  plot  the  kink  angle  for  a  =  0, 
0.4  and  0.8,  where  /?  =  0  and 


£2  =  arctan 


Im  (K*P) 
R  e(Kxhic) ' 


(23) 


The  effect  of  a  on  the  kink  angle  is  stronger  for  small  Q  than  for  large  fi.  In  Fig.  12,  we 
plot  the  kink  angle  for  J?  =  0  and  0.2,  where  a  =  0.4.  Since  for  many  engineering  materials 
is  less  than  0.2  (see  Suga  et  al.,  1988),  the  two  figures  show  that  fS  has  a  weaker  effect  on 
the  kink  angle  than  a.  The  results  which  are  obtained  from  the  linear  variation  of  the 
material  properties  also  support  the  observation. 

The  current  study  provides  a  means  to  assess  crack  deflection  in  a  special  set  of  FGMs, 
simply  nonhomogeneous  materials.  From  the  above  solutions,  it  is  seen  that  material 
gradients  do  have  a  strong  effect  on  the  kink  direction  when  a  crack  is  at  the  central  part 
of  the  FGMs.  The  solutions  may  be  used  in  either  testing  the  fracture  behavior  or  measuring 


3098 


Pei  Gu  and  R.  J.  Asaro 


the  toughness  of  FGMs.  On  the  other  hand,  this  study  can  be  viewed  as  an  initial  effort 
for  crack  deflection  in  real  FGMs :  the  effect  of  microstructure,  i.e.,  the  effect  of  local 
inhomogeneities  around  the  crack  tip  on  both  the  energy  release  rate  and  the  toughness,  is 
up  for  further  investigation. 
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Abstract — This  paper  addresses  failures  near  irregularities  on  the  interface  between  a  film  and  a 
substrate.  Several  boundary  value  problems,  including  two-dimensional  and  three-dimensional 
problems,  involving  inclusions  of  various  shapes  placed  on  the  interface,  are  considered.  The  loading 
is  induced  by  the  lattice  parameter  mismatch  between  the  film  and  substrate.  Stresses  near  the 
interface  and  the  inclusion  boundary  are  of  particular  interest.  The  solutions  show  stress  con¬ 
centration  around  the  inclusion  boundary ;  in  fact,  a  logarithmic  singularity  exists  at  the  intersection 
of  the  inclusion,  film  and  substrate.  Emphasis  is  placed  on  identifying  failures  associated  with  high 
stresses  near  the  inclusion.  A  theoretical  prediction  of  the  misfit  strain  to  cause  adhesion  failure  is 
obtained.  The  driving  force  for  dislocation  emission  from  the  inclusion  is  calculated,  and  it  is  shown 
that  dislocation  emission  from  inclusions  is  favoured  under  a  sufficiently  large  misfit  strain.  ©  1997, 
Elsevier  Science  Ltd.  All  rights  reserved. 


1.  INTRODUCTION 

Thin  films  are  widely  used  in  electronic  components  and  other  advanced  devices.  Because 
high  residual  strains  in  the  films  due  to  either  thermal  expansion  coefficient  mismatch  or 
lattice  parameter  mismatch  between  the  films  and  substrates  can  cause  device  failures 
related  to  dislocation  nucleation  and  interfacial  cracking,  the  study  of  these  mechanisms  is 
a  subject  of  interest.  Many  works  addressing  such  problems  can  be  found  in  the  open 
literature.  Examples  of  recent  development  include  solutions  to  a  large  number  of  interfacial 
cracking  problems  by  Hutchinson  and  Suo  (1991),  a  detailed  treatment  for  threading 
dislocations  in  the  thin  films  by  Freund  (1993),  and  solutions  to  various  stress-related 
problems  in  silicon  technology,  including  film-edge  induced  stresses  and  dislocation  gen¬ 
eration  in  substrate  due  to  these  stresses,  by  Hu  (1991). 

The  interfacial  microstructure  of  a  film/substrate  system  controls  the  performance  of 
the  system  in  more  ways  than  one.  It  appears  that  this  aspect  has  been  overlooked  by  many. 
Recent  experiments  (Perovic  et  al ,  1989)  have  shown  that  jS-SiC  precipitates  are  formed 
on  the  interface  between  the  GevSi]  _x  layer  and  Si  substrate  by  molecular  beam  epitaxy, 
and  dislocations  emit  from  these  precipitates  when  a  critical  particle  size  is  reached.  Also, 
it  is  possible  that  the  surface  of  the  substrate  is  not  a  perfect  plane,  i.e.,  there  are  irregularities 
on  the  substrate.  Inclusions  on  the  interface  are  sources  of  stress  concentration,  and 
therefore  are  the  possible  sites  for  failure  initiation.  Besides  dislocation  nucleation,  adhesion 
failure  occurs  if  stresses  on  the  interface  exceed  its  adhesion  strength.  In  this  paper,  we 
examine  the  stress  concentration  near  these  inclusions,  and  attempt  to  identify  failures 
resulting  from  these  high  stresses.  For  this  purpose,  an  inclusion  on  the  interface  between 
a  thin  film  and  substrate  is  considered,  and  the  loading  arises  from  the  lattice  parameter 
mismatch  between  the  film  and  substrate.  In  this  study,  we  take  the  inclusion,  film  and 
substrate  as  elastically  deforming  bodies,  i.e.,  the  problem  is  solved  by  elasticity  theory. 
The  inclusion  can  be  regarded  as  a  local  inhomogeneity.  So,  it  is  assumed  that  both  the 
thickness  of  the  film  and  the  substrate  are  much  larger  than  the  size  of  the  inclusion.  Under 
these  assumptions,  we  are  able  to  obtain  several  analytical  solutions  for  two-dimensional 
and  three-dimensional  inclusion  problems. 
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The  plan  of  this  paper  is  as  follows.  In  Section  2,  the  two-dimensional  problem,  an 
elliptical  inclusion  on  the  interface  between  a  film  and  substrate,  is  analyzed.  This  idealized 
case  admits  an  analytical  solution  with  a  relatively  simple  form  which  is  obtained  by 
complex  variable  method.  In  Section  3,  the  stresses  on  the  interface  between  the  film  and 
substrate  for  the  two-dimensional  problem  are  examined.  It  is  shown  that  high  stresses  are 
generated  near  the  inclusion ;  in  fact,  the  shear  stress  has  a  logarithmic  singularity  at  the 
intersection  of  the  inclusion,  film  and  substrate.  A  theoretical  prediction  for  the  misfit  strain 
to  cause  adhesion  failure  is  obtained.  Section  4  addresses  the  bending  effect  due  to  the  misfit 
strain.  The  three-dimensional  problem  is  solved  in  Section  5.  Section  6  discusses  dislocation 
emission  from  the  inclusion,  where  the  driving  force  for  dislocation  emission  is  derived, 
and  it  is  shown  that  dislocation  emission  is  favored  for  a  sufficiently  large  misfit  strain. 


2.  TWO-DIMENSIONAL  PROBLEM:  AN  ELLIPTICAL  INCLUSION 

The  mathematical  problem  shown  in  Fig.  1  is  a  film  of  thickness  h  bonded  to  the 
surface  of  a  substrate.  At  the  origin,  there  is  an  elliptical  inclusion  on  the  interface  between 
the  film  and  substrate.  The  major  and  minor  axes  are  a  and  b,  respectively.  The  film  has  a 
lattice  parameter  afi  which  differs  from  that  of  the  substrate  and  inclusion,  a,.  The  inclusion 
and  the  substrate  could  have  different  lattice  parameters,  but  this  possibility  is  not  con¬ 
sidered  here.  The  interface  between  film  and  substrate  and  the  interface  between  film  and 
inclusion  are  assumed  to  be  coherent,  so  that  a  stress  field  is  induced  in  the  film,  substrate 
and  inclusion.  Far  away  from  the  inclusion,  its  solution  is  known:  the  substrate  is  stress 
free  and  the  film  is  subjected  to  biaxial  uniform  tension  or  compression 
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Both  the  film  and  substrate  have  shear  modulus  and  Poisson’s  ratio  vb  while  these  elastic 
constants  for  the  inclusion  are  and  v2.  For  the  semiconductor  layered  systems,  the 
difference  in  elastic  properties  between  the  film  and  substrate  is  usually  small ;  it  is  neglected 
here.  As  stated  in  the  introduction,  the  film  thickness  is  much  larger  than  the  inclusion  size. 

To  solve  the  boundary  value  problem,  one  must  match  lattice  parameters  across  the 
interfaces  between  the  film  and  substrate  and  between  the  film  and  inclusion.  On  the 
interface  between  the  film  and  substrate,  there  is  a  lattice  parameter  mismatch  in  the  x  and 
z  directions  (the  z  axis  is  perpendicular  to  the  x  and  y  axes  shown  in  Fig.  1) ;  but  on  the 
interface  between  the  film  and  inclusion,  there  is  a  lattice  parameter  mismatch  in  all  three 
directions,  x,  y  and  z.  The  solution  to  the  three-dimensional  problem  can  be  obtained  by 
considering  two  two-dimensional  problems  using  an  Eshelby-type  superposition,  as  shown 


Fig.  1.  A  elliptical  inclusion  on  the  interface  between  a  film  and  substrate. 
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(a) 


°V^%°X 


(b) 


Fig.  2.  A  superposition  scheme  for  solving  lattice  parameter  mismatch  problem. 


in  Fig.  2.  First,  we  match  the  lattice  parameters  on  both  interfaces  in  the  x  and  z  directions, 
stretching  the  film  by  e0  in  the  x  and  z  directions.  The  state  of  stress  induced  by  this  process 
is  a  uniform  biaxial  stress,  <jxx  =  oz.  =  <r0.  Next,  a  distribution  of  traction  is  applied  to  the 
interface  between  the  film  and  inclusion  so  as  to  free  it  from  stress,  Fig.  2(b).  In  this  cut 
and  paste  process,  we  have  neglected  the  effect  of  Poisson’s  ratio.  The  stretch  of  the  film 
induced  a  strain  eyy  =  —  2vj/(l  —  v,)£0;  this  strain  reduces  the  lattice  parameter  in  the  y 
direction  of  the  film.  We  have  also  neglected  the  difference  in  the  lattice  parameters  in  the 
y  direction.  These  effects  may  be  neglected  if  both  the  film  and  the  substrate  have  the  same 
lattice  parameter  in  the  y  direction  and  the  Poisson’s  ratio  of  the  film  is  very  small.  We  first 
solve  the  problem  shown  in  Fig.  2(b),  then  the  problem  for  the  mismatch  in  the  y  direction. 

The  problem  shown  in  Fig.  2(b)  may  be  solved  using  the  complex  variable  method 
(Muskhelishvili,  1953).  Stress  and  displacement  fields  are  generated  by  complex  potentials 
<t>(£)  and  ip(£)  such  that : 


c.y.y  +  =  2[$Xc) +  <£(£)]> 

oyy-cxx  +  2\oxy  =  2  [(4>"(c)  +  nt)], 

2fi(u  +  iv)  =  k4>(Q-WW-W)>  (3) 

where  c  =  x+iy;  k  =  3— 4v  for  plane  strain;  k  =  (3  —  v)/(l  +  v)  for  plane  stress.  The 
resultant  force  along  the  arc  from  £  =  £x  to  £  =  c2  can  be  expressed  as 

fx+ifx  =  -im)+mo+W)}\ |:|?.  (4) 

The  resultant  force  due  to  the  traction  acting  on  the  inclusion  boundary  between  £  =  a 
and  c  =  (f0  in  Fig.  2(b)  is 
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Using  a  conformal  mapping,  £  —  R([  +  m/ 0  with  R  =  (a  +  b)/ 2  and  m  =  {a  —  b)j{a-\-b),  the 
elliptical  inclusion  is  mapped  onto  a  unit  circle.  Then,  continuity  conditions  on  the  interface 
between  the  inclusion  and  matrix  (the  film  and  substrate)  can  be  written  as 
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=  0, 

=  <7(t),  (6) 


where 


Imr  ^  0 
Imt  <  0, 


C  = 


/bO+V|) 
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e0/?(l -m). 


(7) 


and  |r|  =  1.  Here,  <£,(£)  and  iAi(Q  are  complex  potentials  defined  for  the  matrix,  and  4>2(C) 
and  tp2{Q  are  complex  potentials  defined  for  the  inclusion. 

The  solution  to  the  problem  is  the  four  potentials  which  satisfy  (6).  Also,  <f> ,(Q  and 
tAi(Q  have  to  give  a  stress  free  state  at  2  =  oo.  The  detailed  procedure  to  solve  the  problem 
is  given  in  the  Appendix,  and  the  solution  to  the  elliptical  inclusion  is  obtained  in  terms  of 
series.  For  a  circular  inclusion,  the  solution  is  a  simple  finite-term  solution,  which  is  given 
by 
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where 


A  = 


ro  +  Kj) 


(T +  k2)(2T +k2  — 1)  2  ’ 


(9) 


and  T  =  ix2lp j.  We  shall  not  write  here  the  expressions  for  02(O  and  \ j/2(Q  since  the  stress 
fields  in  the  film  and  substrate  are  our  interest. 

The  solution  for  the  elliptical  inclusion  can  be  obtained  directly  in  following  special 
cases:  T  —  0,  T  =  1,  and  T  —  oc.  For  T=o o  (rigid  inclusion),  the  solution  is 
0,(0  =  02(O  =  0,(0  =  02( 0  =  0-  When  T  =  3,  the  problem  corresponds  to  a  traction 
(Tody/ds  (d5  is  the  differential  arc  length  of  the  boundary  of  the  ellipse)  acting  on  the  upper 
boundary  of  the  ellipse  in  a  homogeneous  body.  Its  solution  can  be  readily  obtained  by 
integrating  a  point  force  solution  over  the  upper  boundary  of  the  inclusion.  Similarly,  for 
T  =  0.  the  solution  can  be  obtained  by  integrating  a  point  force  acting  on  the  upper 
boundary  of  the  elliptical  hole. 

A  step  irregularity  on  the  substrate  surface  shown  in  Fig.  3  gives  rise  to  a  similar 
problem.  As  before,  both  the  film  and  substrate  have  the  same  elastic  properties;  there  is 
no  lattice  parameter  mismatch  in  the  y  direction,  and  the  Poisson’s  ratio  of  the  film  is  very 
small.  The  solution  to  the  problem  is  obtained  by  the  integration  of  a  point  force  solution 
over  0  <  y  <  a  with  the  point  force  being  <70dy.  The  two  potentials  are  given  by 


<Mc)  =  2rci  J°+  —j  - ia)  lQg(4  -  ia)  ~  C  log  <?], 

<Mc)  =  ~  TZnTTZ  J(^-ifl)log(-:-ifl)-4logc3+  °°  <f[log(<f-ifl)~  logfl.  (10) 
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Fig.  3.  A  step  irregularity. 


We  next  solve  the  problem  for  the  mismatch  in  the  y  direction.  We  view  the  problem 
as  the  one  similar  to  Eshelby’s  internal  problem,  by  assuming  that  there  is  a  uniform 
eigenstrain  in  the  upper  half  of  the  inclusion  (the  portion  above  the  interface  between  the 
film  and  substrate)  s*.  due  to  the  mismatch  in  the  y  direction.  The  eigenstrain  is 


a*  =  s0  1+2 


The  term  involving  v,  represents  the  effect  of  the  Poisson’s  ratio  v,  in  the  cut  and  paste 
process  discussed  before,  and  the  rest  represents  the  effect  of  the  difference  in  lattice 
parameters  in  the  y  direction.  If  there  is  no  lattice  parameter  mismatch  in  the  y  direction, 


£%  =  2-—:- s0. 
"  1—  v, 


The  mathematical  problem  can  be  formulated  by  the  complex  variable  method,  see  List 
and  Silberstein  (1966).  The  continuity  condition  on  the  interface  between  the  inclusion  and 
the  matrix  is  the  same  as  (6),  except  in  this  case  the  right  side  of  the  first  equation  is  g{i) 
and  that  of  the  second  equation  is  zero.  Here, 
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This  shows  that  the  problem  can  be  solved  in  the  similar  way  as  that  in  the  Appendix.  The 
two  potentials  for  a  circular  inclusion  are  given  by 


0,(0  = 


r  H\&%R 

Hm  +  I  2m 


„  in  0-1,  C+I 

—  2+  —  H - ;; — log- — - 

L  L  i— 1 


0,(0  = 


r  /J,e *.R 

T  +  Kr2  2m 


in  0  —  1  C+ 1 

-2+-  +  — bg  — 


—  0'l(Ov  +Ai~, 


where 


A\  —  — 


rq+K,)  Mis*r 

(r  +  K2)(2T+K',-l)  2  ' 


Finally,  we  mention  that  the  analysis  for  the  mismatch  in  the  v  direction  is  the  first  order 
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approximation  since  we  assume  a  uniform  eigenstrain  in  the  inclusion ;  in  the  real  case  the 
eigenstrain  appears  to  be  more  complicated. 


3.  STRESSES  ON  THE  INTERFACE  BETWEEN  THE  FILM  AND  SUBSTRATE 

Having  the  solution  in  Section  2,  we  examine  the  stresses  on  the  interface  between  the 
film  and  substrate  near  an  inclusion.  For  the  most  part  of  this  section,  we  assume  that  the 
film  and  the  substrate  have  the  same  lattice  parameter  in  the  y  direction  and  the  Poisson’s 
ratio  of  the  film  is  small,  so  that  the  solution  for  the  problem  shown  in  Fig.  2(b)  is  a  valid 
solution.  For  a  circular  inclusion,  the  shear  stress  on  the  interface  is 


MP,0) 


<Tq  1 

4  n  1  +rjc, 


3-V 

p4 


log 


P  +  [\ 

p-y 


|  gp  K2 
4n  T + k2 


1+p2,  p  +  l\ 


(16) 


where  p  =  x/a,  and  is  either  greater  than  1  or  less  than  - 1 .  It  has  a  logarithmic  singularity 
at  the  intersection  of  the  film,  substrate  and  inclusion,  p  =  ±  1 .  The  singular  term  is 
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as  p  -*  ±  1.  The  expression  before  the  logarithmic  function  in  the  above  expression  is  a 
logarithmic  stress  intensity  factor,  which  varies  with  T,  k,  and  k2.  For  fixed  ?c,  =  k2  =  2.6 
(plane  strain  with  v,  =  v2  =  0.1),  the  logarithmic  stress  intensity  factor  is  the  largest  when 
T  =  1,  and  it  vanishes  when  T  =  0  (hole)  and  T  =  oo  (rigid  inclusion).  The  singularity  is 
so  weak  that  its  zone  of  dominance  is  very  small.  This  can  be  seen  from  Fig.  4,  where  (16) 
is  plotted  for  several  T’s.  In  this  figure,  the  shear  stress  for  T  =  0  is  larger  than  the  shear 
stress  for  T  =  1  except  the  very  small  portion  near  p  =  1.  Also,  the  shear  stress  decreases 
quickly  to  zero,  and  its  effective  zone  is  about  3-4  times  the  size  of  the  inclusion. 

To  examine  the  geometrical  influence  of  the  inclusion,  the  shear  stress  caused  by  an 
elliptical  inclusion  is  plotted  in  Fig.  5  for  several  values  of  b/a,  where  we  choose  T  =  1  and 
Ki  =  k2  —  2.6  (plane  strain  with  V]  =  v2  =  0.1).  It  shows  that  <jxv  increases  as  b/a  increases. 
This  is  because  the  mismatch  increases  as  b/a  increases.  Also,  the  shear  stress  decays  more 
quickly  for  smaller  b/a. 

An  inclusion  on  the  interface  between  the  film  and  substrate  can  generate  high  shear 
stress.  When  the  inclusion  and  matrix  have  the  same  elastic  properties,  T  =  1  and 
Ki  =  k2  =  2-6  (plane  strain  with  v,  =  v2  =  0.1),  the  shear  stress  at  a  distance  0.1a  ahead  of 


x/a 


Fig.  4.  Shear  stress  on  the  interface  between  the  film  and  substrate  near  a  circular  inclusion. 
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Fig.  5.  The  effect  of  inclusion  geometry  on  shear  stress. 


the  inclusion  (the  shear  stress  is  singular  at  the  intersection  of  the  film,  substrate  and 
inclusion)  is 


<jxy  =  0.23ao  (18) 

for  a  circular  inclusion.  Consider  growing  a  GaAs  film  on  Si  substrate,  and  neglect  the 
difference  in  elastic  properties  between  GaAs  and  Si.  For  these  materials,  e0  —  —0.04  and 
^GaAs  =  85  GPa.  Then,  (18)  gives  axy  =  —0.9  GPa. 

Near  a  circular  inclusion,  the  normal  stress  on  the  interface  between  the  film  and 
substrate  is 


(  m  _  G°  Kl — 1  1  3g°  1  1 

OyyiP’V)  4  2r  +  jc2-l  4  H-Hc,/’ 


(19) 


where  p  =  xja,  and  is  either  greater  than  1  or  less  than  —  1.  We  note  that  it  does  not  have 
a  singularity  as  does  cxy.  For  T  =  1  and  kx  =  k2  —  2.6  (plane  strain  with  Vj  —  v2  =  0.1),  the 
maximum  normal  stress  which  occurs  at  p  =  ±  1  is 


<Tyy  =  —  0.  lcr0.  (20) 

If  Sq  <  0,  the  normal  stress  is  a  peeling  stress.  The  above  expression  shows  that  a  large 
compressive  residual  strain  gives  a  large  peeling  stress  on  the  interface,  which  is  a  driving 
force  for  adhesion  failure  of  the  interface.  If  e0  >  0,  it  is  a  compressive  stress.  For  the 
GaAs/Si  example,  the  maximum  peeling  stress  is  0.4  GPa.  Figure  6  shows  <jyy  vs  p  for 
several  values  of  T. 


Fig.  6.  Normal  stress  on  the  interface  between  the  film  and  substrate  near  a  circular  inclusion. 
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Fig.  7.  Comparison  of  shear  stresses  near  a  circular  inclusion  and  near  a  step  irregularity,  on  the 
interface  between  the  film  and  substrate. 


Near  a  step  irregularity,  the  shear  stress  on  the  interface  between  the  film  and  substrate 
is 
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where  p  =  xja,  and  p  >  0.  The  shear  stresses  on  the  interface  near  a  step  irregularity  and 
near  a  circular  inclusion  are  compared  in  Fig.  7.  As  expected,  it  shows  that  the  step 
irregularity  produces  larger  shear  stress  on  the  interface. 

If  there  are  no  inclusions  and  step  irregularities  on  the  interface,  there  are  no  stresses 
acting  on  the  interface.  Therefore,  adhesion  failure  may  not  occur.  If  the  residual  strain  s0 
in  the  thin  film  causes  adhesion  failure  due  to  inclusions  or  step  irregularities,  we  estimate 
it  from  above  solutions.  Suppose  that  the  interfacial  adhesion  strength,  i.e.,  the  summation 
of  all  atomic  bonding  forces  between  the  adjoining  material  surfaces  for  an  ideal  planar 
interface,  is  known,  and  let  the  critical  normal  stress  be  ab  and  the  critical  shear  stress  be 
Tb.  Then,  from  (18)  and  (20),  the  mismatch  stress  cr0  to  cause  adhesion  failure  is  given  by 
the  smaller  of 


I*oI=4.35t„ 

(22) 

for  compressive  <r0.  For  tensile  cr0,  this  prediction  is 

<70  =  4.35t*.  (23) 

In  writing  (22)  and  (23),  we  have  assumed  that  adhesion  failure  caused  by  shear  occurs  at 
a  distance  0.1  a  ahead  of  the  inclusion.  Besides  calculating  atomic  bonding  forces,  the 
interface  adhesion  strength  can  be  measured  by  a  number  of  methods  which  have  been 
discussed  in  Alexopoulos  and  O'Sullivan  (1990).  Due  to  the  difficulty  in  handling  these 
tests,  the  data  obtained  by  these  methods  are  usually  considered  to  be  qualitative  in  nature. 

It  is  seen  that  either  a  large  residual  strain  or  poor  interface  bonding  could  lead  to 
adhesion  failure.  In  the  GaAs/Si  example  discussed  above,  the  film  is  subjected  to  com¬ 
pressive  residual  stress  and  the  normal  stress  on  the  interface  is  a  peeling  stress.  If  we 
assume  that  the  failure  mode  is  the  adhesion  failure  by  the  peeling  stress,  the  criterion 
predicts  that  the  adhesion  failure  occurs  if  ob  is  less  than  0.4  GPa.  The  development  of  the 
initial  failure  could  lead  to  large  scale  fracture,  such  as  film  buckling.  The  mechanics  of 
film  buckling  has  been  studied  by  Evans  and  Hutchinson  (1984),  and  the  cracking  and 
buckling  processes  of  a  film/substrate  system  have  been  characterized  in  terms  of  the 
behavior  of  the  film  and  substrate,  and  the  interface  bonding  by  Evans  el  al.  (1988). 
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The  effect  of  the  mismatch  in  the  y  direction  can  be  estimated  using  the  solution  for  a 
circular  inclusion  given  by  (14)  and  (15).  The  shear  stress  on  the  interface  between  the  film 
and  substrate  near  the  inclusion  is 


O. xv(p,  0) 


/fi£,n 


271  1+FK] 


6  3  —  p2  p-Fl 

-  -  +  — r-iog  — 


p~K 
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2  1  +p2  p  + 1 

-“  +  — “log 
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P~  1 


2jl  T  +  K; 

where  p  =  xja,  and  is  either  greater  than  1  or  less  than  —  1 ;  the  normal  stress  is 
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ffvr(P»0)  =  P\£ 


_ 1  ip^tr  r  1 

2T + k'2  —  1  p2  2  1+Hc,  p*' 


(24) 


(25) 


We  examine  the  effect  of  the  Poisson’s  ratio  using  (12)  for  £*..  Since  the  stresses  are 
proportional  to  the  Poisson’s  ratio  v,,  they  can  be  neglected  if  the  Poisson’s  ratio  is  small. 
The  larger  the  Poisson’s  ratio,  the  larger  its  effect  on  the  stress  field.  Few  numerical  examples 
for  the  stresses  are  given  below  when  v,  is  small,  by  taking  T  =  1  and  tc,  =  k2.  At  p  =  1, 
the  normal  stress  is  0.006  <x0  for  V)  =  0.01,  and  0.06  a o  f°r  vi  =  0.1.  At  p  =  1.1,  the  shear 
stress  is  0.001  o0  for  v,  =  0.01,  and  0.01  <r0  for  v,  =  0.1.  It  appears  that  the  normal  stress  is 
larger  than  the  shear  stress  for  the  mismatch  in  the  y  direction. 


4.  BENDING  EFFECT 

If  the  film  and  substrate  have  finite  thicknesses,  the  solution  described  in  the  preceding 
sections  must  be  corrected.  Four  self-equilibrating  bending  moments  are  applied  to  the  film 
and  substrate,  and  two  of  them  in  the  x-y  plane  are  shown  in  Fig.  8.  We  assume  the 
thickness  of  the  substrate  hs  is  much  larger  than  that  of  the  film  h,  and  the  thickness  of  the 
film  is  much  larger  than  the  size  of  the  inclusion  R=(a+b)/ 2.  The  magnitude  of  these 
bending  moments  depends  on  the  thicknesses  of  film  and  substrate,  and  is  given  by 

M  =  \  hhsaQ .  (26) 

These  bending  moments  act  at  a  distance  hs/2  from  the  bottom  of  the  substrate.  For  T  =  1 
and  k,  =  k-2,  M  does  not  cause  shear  and  normal  stresses  at  the  interface  between  the  film 
and  substrate,  and  the  stress  and  strain  fields  in  the  film  and  substrate  are  given  by  beam 
theory.  For  T  /  ] ,  the  bending  moments  induce  stresses  at  the  interface,  which  we  calculate 
as  follows. 

Note  that  the  two  moments  in  the  direction  do  not  cause  stresses  on  the  interface  if 
the  inclusion  and  the  matrix  have  the  same  curvature  in  the  y-z  plane  (the  -  axis  is 
perpendicular  to  the  x  and  y  axes  shown  in  Fig.  1).  Consequently,  it  is  only  necessary  to 
consider  the  deformation  caused  by  the  two  moments  shown  in  Fig.  8.  If  we  use  the 
displacement  field  of  a  homogeneous  problem  (T  =  1,  k,  =  x-2)  as  a  solution  to  the  inhomo¬ 
geneous  problem  there  is  a  traction  jump  on  the  inclusion  boundary.  The  solution  may  be 


Fig.  8.  Bending  moments  acting  on  the  film  substrate  system. 
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corrected  by  superposing  an  equal  and  opposite  traction  jump  on  the  inclusion  boundary. 
The  traction  jump,  expressed  as  a  resultant  force  on  the  inclusion  boundary,  is 

,  .  .  ,  .  1  „  1 

9{t)  =  1C,T-+1C,  —  —  2lCi  +  c2t-c2-, 

T*  1 

1  M 

Cy  =  g^2(l -/»)-(! -T)  —  , 

c2  =  -Iy?(i -,„)(i -T)^,  (27) 

where  1  =  (hs+h)3/\2.  For  simplicity,  we  have  taken  jc,  =  k2  to  derive  (27).  When  hs  ->  oo, 
ff(r)  =  0,  which  states  the  fact  that  a  substrate  with  infinite  thickness  does  not  induce 
bending  curvature.  The  function  g( t)  plays  the  same  role  as  that  in  Section  2,  so  that  the 
solution  can  be  obtained  by  the  same  way  as  that  in  the  Appendix.  For  a  circular  inclusion, 
the  solution  is 


<M0  = 


rK-+ulc,C2 


.  K  .  1  1  —  K  1  1,  ... 

_  r+K1C|  “  2r+K- 1 C2 ( "  c 


(28) 


where  k  —  k,  =  k2.  Then,  the  normal  and  shear  stresses  on  the  interface  between  the  film 
and  substrate  near  the  inclusion  are 
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(29) 


The  condition,  h,  hs  »  must  be  satisfied  for  the  solution  to  be  valid  as  explained  at  the 
beginning  of  the  section.  The  normal  stress  is  in  compression  for  a  soft  inclusion  (T  <  1), 
and  the  shear  stress  is  proportional  to  a  coefficient  in  (29),  i.e.. 


0) 


3Rh 


(30) 


Suppose  that 


hs  >  10/r  >  1000/?, 


we  have  from  (30)  that 


<rxy(j>,  0) 


9 

80000  <T°‘ 


(31) 


(32) 


The  shear  stress  may  be  neglected. 

5.  THREE-DIMENSIONAL  PROBLEM 

In  general,  the  problem  of  an  inclusion  on  the  interface  between  the  film  and  substrate 
is  likely  to  be  a  three-dimensional  problem.  The  analytical  solution  to  a  three-dimensional 
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Fig.  9.  Three-dimensional  problem. 


inclusion  problem  is  difficult  and  available  solutions  can  be  found  in  Mura  (1987).  In  this 
section,  by  assuming  that  the  elastic  properties  of  the  inclusion,  film  and  substrate  are  the 
same,  we  find  an  exact  solution  to  the  three-dimensional  inclusion  problem  for  the  case 
that  the  loading  is  lattice  parameter  mismatch  between  the  film  and  substrate.  As  before, 
we  assume  that  the  film  thickness  is  much  larger  than  the  size  of  the  inclusion. 

The  three-dimensional  problem  can  be  decomposed  into  two  problems.  One  represents 
the  lattice  parameter  mismatch  in  the  x  and  y  directions,  and  the  other  represents  the 
mismatch  in  the  z  direction  (the  x,  y,  z  axes  are  shown  in  Fig.  9,  and  the  film  thickness 
direction  is  designated  as  z).  As  in  Section  2,  we  shall  assume  that  there  is  no  lattice 
parameter  mismatch  in  the  z  direction  and  the  Poisson’s  ratio  of  the  film  is  very  small.  In 
other  words,  the  mismatch  in  the  z  direction  can  be  neglected.  Following  the  argument 
given  in  Section  2,  the  solution  to  the  three-dimensional  problem  can  be  obtained  by 
integrating  a  point  force  solution  over  the  inclusion  boundary.  To  do  this,  we  write  the 
solution  of  a  point  force  in  the  following  form 

pk 

4(x, y,z)  =  tfjix-XQty-yotZ-Zo),  (33) 

where  indices  Uj  and  k  vary  from  1  to  3,  designating  the  x,  y  and  z  in  turn ;  Pk  is  a  vector 
force  in  the  k  direction.  The  function  <r*v  has  the  dimension  1 /[length]2.  Equation  (33) 
represents  the  stress  field  due  to  a  point  force  in  the  direction  k  acting  at  (x0,  y0,  z0) .  The 
stresses  induced  by  an  inclusion  or  a  step  irregularity  can  be  written  as 


<?ij(x9y9z)  = 


8tt(1  —  v) 


(*oO'0’rO>e5 


nx6}j(x-x0,y-y0,z-29)  dS 


+  ny6l(x  -  x0 ,  j  -  t'o  ,  z  -  z0 )  dS 

J(.Y0.v0.r0)eS 


(34) 


where  S  is  the  inclusion  boundary  that  is  surrounded  by  the  film,  or  the  vertical  part  of  the 
step  irregularity,  and  n  =  (nY,  «v,  n.)  is  the  unit  outward  normal  to  the  inclusion  boundary 
at  (x0,j*o,Zo).  For  the  two  normal  stresses  parallel  to  the  interface  between  the  film  and 
substrate,  axx  and  <rvr,  the  residual  stress  er0  has  to  be  added  to  obtain  the  total  stresses.  The 
two  surface  integrals  in  (34)  can  not  be  evaluated  analytically  for  an  arbitrary  inclusion 
shape.  They  cannot  even  be  evaluated  exactly  for  some  simple  shapes  like  a  sphere.  However, 
the  surface  integrals  can  be  reduced  to  line  integrals  to  simplify  them  in  many  cases.  Two 
examples  will  be  given  below. 

The  first  example  is  a  corner  step  irregularity,  as  shown  in  Fig.  10.  The  corner  step 
irregularity  occupies  the  region,  x  ^  0,  v  ^  0  and  0  ^  z  ^  a.  Note  that,  if  we  extend  the 
above  step  irregularity  region  in  the  y  direction  and  let  it  occupy  x  ^  0,  —  oo  <  y  ^  oo  and 
0  ^  z  ^  af  it  becomes  the  two-dimensional  problem  which  was  solved  in  Section  2.  The 
purpose  of  this  example  is  to  examine  the  effect  of  a  corner  on  the  stress  distribution.  In 
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Fig.  10.  A  corner  step  irregularity. 


this  case,  the  integrations  in  (34)  can  be  evaluated  analytically.  The  two  shear  stresses  on 
the  plane  z  —  0  are  given  by 


axz(x>y,  0)  =  • 


<7o 


a2y 


871(1— V)|_  {xl+al)Jx2+y2+a2  x2+a2 


+  (1  —  2v)  log 
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gp 


a2x 


8ti(1  -v)[  (y2+aVx2-fy2+a2  y2+a2 


+  (1— 2v)log 


■x+y^+y+a2 


(35) 


The  shear  stress  oy.  in  (35)  has  a  logarithmic  singularity  on  the  straight  line,  x  =  0  and 
y  <  0,  as  in  the  two-dimensional  case.  The  corner  introduces  a  stress  cy.,  which  is  zero  in 
the  two-dimensional  case.  Far  away  from  the  corner,  its  influence  vanishes  so  that  <jy. 
approaches  zero  as  y  -*  ±  oo.  A  plot  of  ay.  vs  y  when  x  =  0.01a  is  given  in  Fig.  1 1.  We 
observe  that  the  a,.,  is  negligible  for  |x|  >  5a.  On  the  other  hand,  ax.  approaches  the  two- 
dimensional  solution  given  in  Section  2  as  y  ->  —  oo,  and  vanishes  when  +oo.  The 
variation  of  ax;  along  the  y  axis  when  x  =  0.01a  is  also  shown  in  Fig.  11.  The  two  shear 
stresses  at  y  =  0,  z  =  0  are  given  by 


a.v;(x,  0, 0) 
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Fig.  1 1 .  Shear  stresses  near  a  corner  step  irregularity. 
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+  (]  —  2v)  log 


X  +  y/x*j 

2x  J 


(36) 


for  x  >  0.  Comparing  trx:(x,  0, 0)  with  <rx.(x,  -  oo,  0),  we  see  that  there  is  a  simple  relation¬ 
ship  between  them, 

ax;  (x,  —  oo ,  0)  =  2a  x.  (x,  0, 0) .  (37) 

This  tells  us  that  by  traveling  from  y  =  —  oo  to  y  =  0,  ax:  reduces  to  half  of  its  value  at 
v  =  —  oo,  and  ay.  rises  from  zero  to  the  one  given  in  (36).  Also,  we  have  another  relationship, 

<r2.(x, 0, 0)  +  cry.(x, 0, 0)  <  ffl-(x,  -  oo, 0)  (38) 

for  x  >  0.  This  says  that  on  the  interface  between  the  film  and  substrate  the  total  shear 
stress  at  (x,  0, 0)  is  less  than  the  total  shear  stress  at  (x,  -  oo,  0),  because  the  corner  relieves 
stress  concentration. 

The  second  example  is  a  circular  plate  above  the  interface  of  the  film  and  substrate,  as 
shown  in  Fig.  9.  Both  the  thickness  and  radius  of  the  circular  plate  are  a.  The  axisymmetry 
of  the  problem  renders  <7r:  (r,  8  and  z  are  used  instead  of  x,  y  and  z)  is  the  only  nontrial 
shear  stress,  which  on  the  interface  is  given  by 
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(r2  —  2arcos9+a2y,!-j 


d0. 


(39) 


after  evaluating  the  integrals  in  (34)  along  the  z  direction.  Further  derivation  shows  that 
the  line  integrals  in  above  expression  can  be  written  in  terms  of  elliptical  integrals.  The  arz 
has  a  logarithmic  singularity  at  r  ~  a ,  and  a  plot  of  <rr:  vs  r  is  shown  in  Fig.  12. 

We  can  also  give  the  expressions  for  other  stress  components  in  the  above  two  examples 
by  evaluating  the  integrals  in  (34).  The  inclusion  shape  has  a  strong  effect  on  the 
stresses  near  the  inclusion.  For  the  circular  plate,  <rr.(l.la, 0,0)  =  0.186a0  and  azz 


Fig.  12.  Shear  stress  near  a  circular  plate  inclusion. 


962 


Pei  Gu  ef  a!. 


(a.  0, 0)  =  -  0.02<xo,  whereas  for  the  two-dimensional  circular  inclusion  in  Section  2,  the  two 
stresses  on  the  interface  between  the  film  and  substrate  are  0.23cro  and  —  0.1ct0>  respectively. 


6.  DISLOCATIONS  IN  THIN  FILMS 

Over  the  past  few  decades,  the  so-called  threading  dislocations  which  form  in  epitaxial 
layers  used  for  semiconductor  materials  have  received  attention.  Threading  dislocations 
occur  when  the  thickness  of  an  epitaxial  layer  reaches  a  critical  value.  The  critical  thickness 
has  been  treated  by  a  number  of  authors,  including  Matthews  and  Blakeslee  (1974),  People 
and  Bean  (1985)  and  Freund  (1990).  But  the  mechanism  responsible  for  forming  such 
dislocations  has  not  yet  been  fully  understood.  Eaglesham  et  al.  (1989)  have  discussed  the 
diamond  defect  as  a  source  for  dislocation  nucleation.  Dodson  (1988)  has  argued  that 
the  stress  concentration  around  clusters  of  impurity  atoms  is  a  good  source  to  generate 
dislocations.  The  experiments  by  Perovic  et  al.  (1989)  have  shown  that,  for  films  grown  by 
molecular  beam  epitaxy,  there  are  precipitates  on  the  interface  between  the  film  and 
substrate,  and  V-shaped  threading  dislocations  are  emitted  from  the  heterogeneous  particles 
if  their  size  exceeds  critical  value,  while  particles  below  the  critical  size  are  coherent  with 
surrounding  matrix.  Among  the  two  types  of  heterogeneous  precipitates  on  the  interface, 
it  was  found  that  yS-SiC  precipitates  tend  to  generate  dislocations,  whereas  SiOx  precipitates 
are  dislocation-free  (Hull,  1986)  because  they  are  in  an  amorphous  state  which  eliminates 
misfit  strains.  In  this  section,  we  consider  these  observations  in  the  light  of  our  solutions  to 
the  state  of  stresses  near  inclusions  in  strained  layer  systems. 

We  examine  the  shear  stress,  which  provides  the  glide  force  for  dislocation  motion,  on 
possible  slip  planes  around  the  inclusion.  Two  representatives  of  the  slip  planes  shown  in 
Fig.  13  are  chosen  for  the  two-dimensional  problems  (circular  inclusion  and  step  irregu¬ 
larity)  whose  solutions  have  been  obtained  in  Section  2.  The  orientation  of  these  slip  planes 
is  chosen  here  to  form  the  angle,  0  =  sin-1  (1/^/3),  with  the  interface  between  the  film  and 
substrate.  We  neglect  the  difference  in  elastic  properties  among  the  inclusion,  film  and 
substrate.  Figure  14  shows  the  shear  stress  for  the  circular  inclusion,  where  r/a  is  the 
distance  from  the  inclusion ;  it  has  the  same  trend  for  the  step  irregularity.  The  shear  stress 
at  the  inclusion  boundary  on  slip  plane  #1  for  the  step  irregularity  has  a  logarithmic 
singularity,  and  for  the  circular  inclusion  it  is  nearly  twice  as  large  as  that  far  away  from 
the  inclusion,  the  homogeneous  solution.  The  high  shear  stress  state  provides  a  good 
condition  to  generate  dislocations.  The  shear  stress  on  slip  plane  #2  is  somewhat  different 
from  slip  plane  #1.  Although  there  is  a  logarithmic  singularity  at  the  inclusion  boundary, 
the  shear  stress  changes  sign  at  a  distance  very  close  to  the  inclusion,  and  increases  to  the 
homogeneous  solution  as  the  distance  increases.  The  shear  stress  near  the  inclusion  on  any 


Slip  plane  #1 


Fig.  13.  Two  possible  slip  planes  which  intersect  the  circular  inclusion  and  the  step  irregularity. 
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slip  plane,  which  is  between  slip  planes  #  1  and  #  2  and  is  parallel  to  these  two  slip  planes, 
is  less  than  that  on  slip  plane  #  1 . 

If  the  inclusion  is  admitted  as  a  source  for  dislocation  nucleation  in  a  thin  film,  the 
driving  force  for  the  emitted  dislocation  can  be  derived  as  follows.  The  energy  dissipation 
associated  with  the  dislocation  motion  is  expressed  as 


W{S)=-]- 


afin fate¬ 
's 


afjtijb i  ds, 

5 


(40) 


where  S  is  the  glide  plane ;  tj,  is  the  outward  normal  of  the  left  portion  of  the  film  divided 
by  the  slip  plan;  the  Burgers  vector  6,  is  the  displacement  of  the  left  side  of  the  glide  plane 
minus  that  of  the  right  side  of  the  glide  plan.  In  the  expression,  the  stress  field  afj  is  due  to 
the  dislocation,  and  a*  is  due  to  the  applied  load,  namely  the  lattice  parameter  mismatch. 
Define  the  driving  force  G  from  the  time  rate  of  the  energy  dissipation  as 

H/(S)  =  Gv.  (41) 

where  v  the  speed  of  the  emitted  dislocation.  The  driving  force  is  derived  in  the  same  way 
as  that  in  Freund  (1990)  for  threading  dislocation,  and  is  given  by 

G=  -\o?M)nfa-o?j{B)nfa,  (42) 

where  a f(A)  and  o?,{B)  denote  the  corresponding  stresses  at  the  inclusion  boundary  and 
the  position  of  the  emitted  dislocation,  respectively.  Note  that  since  the  difference  in  elastic 
properties  between  the  inclusion  and  the  matrix  is  neglected,  there  is  no  image  force  due  to 
the  mismatch  of  elastic  properties  acting  on  the  dislocation.  The  first  term  in  (42)  can  be 
calculated  from  the  solution  of  a  dislocation  in  a  infinite  body,  and  the  second  term  can  be 
calculated  from  the  solution  obtained  in  Section  2.  We  finally  obtain 
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where  r  is  the  distance  between  the  dislocation  and  the  inclusion  boundary,  a  is  the  size  of 
the  inclusion,  b  is  the  magnitude  of  the  Burgers  vector,  and  T(r/a)  is  a  nondimensional 
function  obtained  from  the  stresses  <^(5).  Given  the  inclusion  size,  the  magnitude  of  the 
Burgers  vector  and  the  misfit  strain,  the  distance  from  the  inclusion  at  which  G  =  0 
(equilibrium  point),  namely  r*,  is  determined  from  (43).  For  r  >  r*,  G  >  0;  for  r  <  r*, 
G  <  0.  If  the  dislocation  is  formed  at  the  distance  less  than  r*  from  the  inclusion,  it  will  be 
driven  back  to  the  inclusion  boundary ;  otherwise,  it  will  be  emitted.  We  take  the  circular 
inclusion  as  an  example,  and  consider  slip  plane  #1.  For  b  —  0.4  nm  and  e0  =  0.04, 
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r*  =  0.458  nm  if  a  ~  100  nm ;  r*  =  0.473  nm  if  a  =  10  nm ;  r*  =  0.673  nm  if  a  =  1  nm. 
When  £0  is  increased  to  0.06,  r*  is  0.305  nm,  0.31 1  nm  and  0.415  nm  for  the  above  three 
cases,  respectively.  The  larger  the  misfit  strain  and  the  inclusion  size,  the  smaller  the  critical 
distance  r*.  When  r*  is  less  than  the  core  cut-off  radius,  we  have  the  case  that  the  dislocation 
formed  at  the  inclusion  boundary  is  emitted  spontaneously  (see  Rice  and  Thomson  (1974) 
for  dislocation  emission  from  a  crack  tip).  The  numerical  results  here  show  that  dislocation 
emission  is  favored  for  large  misfit  strain.  For  inclusions  in  shapes  other  than  the  two  we 
have  treated  here,  the  driving  force  may  be  obtained  in  the  same  manner.  Zhang  and  Yang 
(1994)  have  studied  the  process  of  dislocation  emission  from  a  spherical  cavity  on  the 
interface  between  the  film  and  substrate. 


7.  CONCLUDING  REMARKS 

Solutions  are  presented  for  several  boundary  value  problems  involving  various  shaped 
inclusions  located  on  the  interface  between  a  thin  film  and  a  substrate.  The  film  and  the 
substrate  have  different  lattice  parameters.  The  solutions  exhibit  stress  concentrations 
around  the  inclusions.  There  is  a  logarithmic  stress  singularity  at  the  intersection  of  the 
film,  substrate  and  inclusion.  Moreover,  the  shear  stress  around  the  inclusion  can  be  nearly 
twice  as  large  as  the  remote  field  given  by  the  homogeneous  solution.  High  stresses  around 
inclusions  and  on  the  interface  between  film  and  substrate  could  cause  dislocation  nucleation 
and  adhesion  failure.  A  theoretical  prediction  of  the  misfit  strain  to  cause  adhesion  failure 
is  obtained.  The  driving  force  for  dislocation  emission  from  the  inclusion  is  calculated,  and 
it  is  shown  that  dislocation  emission  from  inclusions  is  favored  for  sufficiently  large  misfit 
strain.  The  results  suggest  that  small  inclusions  (compared  to  film  thickness)  on  the  interface 
between  a  film  and  a  substrate  can  be  the  sources  for  failure  initiation  in  the  layered  system. 
These  failures  can  be  precursors  to  failures  on  a  larger  scale :  adhesion  failure  can  lead  to 
film  buckling;  and  dislocation  emission  eventually  causes  threading  dislocations  as  shown 
by  Zhang  and  Yang  (1994)  for  a  spherical  cavity. 
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In  obtaining  (A.9)  and  (A. 10),  we  have  made  use  of  (A.4).  Also,  we  have  replaced  the  variable  f  by  in  (A.10). 
By  (A.9) -(A.  10),  an  ordinary  differential  equation  for  0;:( f)  is  obtained, 
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Using  (7),  we  evaluate  the  integrals  in  (A. 6)  and  (A. 8)  to  obtain 
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Generally,  the  solution  for  </>22(£)  in  (A.l  1)  is  expressed  in  terms  of  series.  We  write 
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where  /„  (n  =  1.2, . . .)  are  the  known  coefficients  from  the  series  expansion  of  (A.13).  We  obtain  <f>22(C)  by 
substituting  (A.  14)  into  (A .11)  and  comparing  the  coefficients  of  both  sides. 

By  (A. 5. 1 )  +  (A. 5.2)  x  2/i2,  the  potential  <£,(£)  is  obtained  as 
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Evaluating  the  integrals  in  (A.  15)  and  (A.  16),  we  have  the  final  result, 
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The  two  potentials  for  the  inclusion,  <p2{ 0  and  ^2(e),  can  also  be  obtained. 
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Abstract — A  semi-infinite  crack  in  a  strip  of  an  isotropic,  functionally  graded  material  under  edge 
loading  and  in-plane  deformation  conditions  is  analyzed.  Mixed  mode  stress  intensity  factors  are 
analytically  solved  for  up  to  a  numerically  determined  parameter.  The  effects  of  material  gradients 
on  the  mode  I  and  mode  II  stress  intensity  factors  and  the  phase  angle  used  to  measure  mode  mixity 
are  determined.  The  solution  is  extended  to  the  case  where  the  strip  is  made  of  an  orthotropic, 
functionally  graded  material.  These  results  are  applied  to  solve  a  four-point  bending  specimen 
configuration  that  may  be  used  to  test  the  fracture  behavior  of  functionally  graded  materials.  The 
nature  of  the  crack  tip  fields  and  possible  fracture  criterion  for  functionally  graded  materials  are 
discussed.  Copyright  ©  1996  Elsevier  Science  Ltd 


1.  INTRODUCTION 

Functionally  graded  materials  (FGMs)  to  be  used  as.  inter  alia ,  superheat-resistive  materials 
have  promised  attractive  applications  in  furnace  liners,  space  structures,  and  fusion  reactors. 
FGMs  consist  of  two  distinct  material  phases,  such  as  ceramic  and  metal  alloy  phases,  and 
are  the  mixture  of  them  such  that  the  composition  of  each  changes  continuously  along  one 
direction.  The  change  in  microstructure  induces  chemical,  material,  and  microstructural 
gradients,  and  makes  functionally  graded  materials  different  in  behavior  from  homogeneous 
materials  and  traditional  composite  materials  (Yamanouchi  et  al ,  1990;  Holt  et  al ,  1993). 
These  materials  are  tailorable  in  their  properties  via  the  design  of  the  gradients  in  chemistry 
and  microstructure  that  is  possible  within  them. 

Experiments  have  shown  that  cracks  occur  in  functionally  graded  materials  (see  above 
references)  although  the  absence  of  sharp  interfaces  does  alleviate  problems  with  interface 
fracture.  For  cracks  in  this  type  of  material,  stress  intensity  factors  are  affected  by  the 
material  gradients.  Moreover,  the  fracture  modes  of  the  cracks  in  FGMs  are  inherently 
mixed  when  they  are  not  parallel  to  the  direction  of  material  property  variation,  i.e .,  there 
are  typically  both  normal  and  shear  tractions  ahead  of  the  crack  tips  because  of  the  non¬ 
symmetry  in  the  material  properties.  To  characterize  the  material,  fracture  toughness  data 
is  required.  To  obtain  the  fracture  toughness  data,  stress  intensity  factors  for  specimens 
subjected  to  variable  external  loads  are  needed.  Most  previous  works  on  FGM  crack 
configurations  have  concentrated  on  finite  crack  problems,  e.g.,  Delale  and  Erdogan  (1983, 
1988)  and  Noda  and  Jin  (1993)  have  analyzed  a  finite  crack  in  a  plate  subjected  to 
mechanical  and  thermal  loads.  A  semi-infinite  crack  in  an  interlayer  between  two  dissimilar 
materials  was  considered  by  Yang  and  Shih  (1994),  and  they  obtained  an  approximate 
solution  from  a  known  bimaterial  solution.  We  consider  herein  a  semi-infinite  crack  in  a  strip 
of  an  isotropic,  functionally  graded  material  under  edge  loading  and  in-plane  deformation 
conditions.  Stress  intensity  factors  for  the  crack  tip  are  obtained.  The  solution  is  analytical 
up  to  a  parameter  which  is  obtained  numerically.  The  solution  is  extended  to  the  case  where 
the  strip  is  made  of  an  orthotropic,  functionally  graded  material.  The  results  are  applied  to 
analyze  a  four-point  bending  specimen  configuration  that  may  be  used  to  test  the  fracture 
behavior  of  functionally  graded  materials.  The  mode  III  stress  intensity  factor  in  the  cracked 
plate  subjected  to  anti-plane  deformation  is  obtained.  The  nature  of  the  crack  tip  fields  and 
possible  fracture  criterion  for  functionally  graded  materials  are  discussed. 
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The  main  emphasis  here  is  to  analyze  fracture  behavior  in  materials  that  possess 
continuously  graded  microstructures.  On  the  other  hand,  the  physical  picture  developed 
for  the  cracked  microstructure  actually  provides  a  more  realistic  model  for  cracks  along 
interfaces,  in  general,  at  least  for  those  that  have  any  but  atomic-scale  width. 


2.  FIELD  EQUATIONS  FOR  ISOTROPIC  MEDIA 
In  this  study,  we  take  the  elastic  properties  to  be  of  following  exponential  forms 


E(y)  =  E0e\ 

v'(>’)  =  v0(l+£.v)  ev',  (1) 


where  y  and  £  are  material  constants  representing  the  material  gradients ;  E0  and  v0  are  the 
values  of  these  elastic  properties  at  y  =  0.  For  plane  stress  problems,  E'(y)  -  E(y)  and 
v'(v)  =  v(y),  where  E(y)  and  v(v)  are  Young's  modulus  and  Poisson’s  ratio,  respectively; 
for  plane  strain  problems,  E'(y)  =  E()i}/[\-v(y)2]  and  v'(j')  =  v(v)/[l -v(v)].  The  par¬ 
ameters  y  and  £  have  a  dimension  [length] These  forms  for  the  material  properties  have 
been  previously  used  by  Delale  and  Erdogan  (1988)  and  Noda  and  Jin  (1993) ;  they  provide, 
on  the  one  hand,  analytical  flexibility  and  yet  lead  to  somewhat  simple  forms  for  the  field 
equations.  The  shear  modulus,  pt(y)9  relates  to  Young’s  modulus  and  Poisson’s  ratio  by 


tiy)  = 


goo 
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Using  (1)  and  (2),  the  stress  function  4>(.\%y)  defined  in  the  same  way  as  that  for 
homogeneous  materials,  i.e.,  stresses  are  obtained  from  the  second  derivatives  of  the  stress 
function,  satisfies  the  following  equation 
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For  a  traction  problem,  the  solution  satisfies  (3)  and  boundary  conditions.  The  material 
constant  y  enters  the  stress  field  of  the  traction  problem,  whereas  the  other  material 
parameters  £0,  v0  and  £  do  not.  In  (3),  the  fourth-order  differential  terms  do  not  involve  y. 
and  constitute  the  biharmonic  equation,  which  is  the  equation  for  homogeneous  materials. 
By  dimensional  analysis,  the  stress  field  has  the  following  generic  form 


c„(x,y)  =  TG*[yh.j;jr± 


(4) 


where  /,  j  =  1,2;  Tis  a  representative  stress  magnitude;  h  is  a  characteristic  length  in  the 
problem;  if  is  the  group  of  lengths  which  represents  the  geometry  of  the  problem.  This 
differs  from  the  case  of  a  homogeneous  material  in  which  material  properties  do  not  enter 
the  stress  field  of  a  traction  problem,  and  also  differs  from  the  case  of  a  bimaterial  in  which 
Dundurs’  parameters  (Dundurs,  1969)  measuring  the  material  mismatch  enter  the  stress 
field  of  a  traction  problem. 

The  parameters  yh  in  the  solution  is  dependent  on  the  thickness  of  a  functionally 
graded  material,  L,  the  Young's  moduli  at  the  upper  and  the  lower  boundaries  of  the 
material,  Eu  and  E'h  and  the  characteristic  length,  /?.  From  (1). 
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It  is  seen  that  yh  is  proportional  to  hjL ,  and  increases  logarithmically  with  the  ratio  of  the 
two  Young’s  moduli.  For  example,  if  hjL  =  1,  yh  =  0.35  for  E'JE]  =  0.5,  and  yh  =  0.97 
for  E'JE J  =  7.  The  choice  of  h  is  arbitrary.  If  h}  and  h2  denote  two  choices  for  the  charac¬ 
teristic  length,  the  corresponding  stresses  obtained  satisfy  the  relation 
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3.  CRACK  TIP  FIELDS 

A  brief  review  of  the  crack  tip  fields  in  functionally  graded  materials  is  given  in  this 
section.  Consider  a  crack  in  a  strip  of  a  functionally  graded  material,  as  illustrated  in  Fig. 
1.  Stresses  near  the  crack  tip  have  a  square-root  singularity,  and  singular  terms  of  the 
stresses  (Jin  and  Noda,  1994)  are  of  the  form 


=111 
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where  /,  j=  1,  2;  r  and  6  are  the  polar  coordinates  shown  in  Fig.  1.  The  dimensionless 
angular  functions  a!j(0)  and  cf(6)  are  the  same  as  those  for  homogeneous  materials. 
This  can  be  easily  proved  by  expanding  the  stress  function  as  S/=or'+p®/0,y)>  sub¬ 
stituting  the  series  into  (3).  The  resulting  equation  for  <£o(0,  y)  and  the  eigen-value  problem 
used  to  determine  p  do  not  involve  y  and  are  the  same  as  those  for  homogeneous  materials. 
In  fact,  for  any  form  of  material  properties  and  any  orientation  of  the  crack,  the  highest 
order  differential  terms  in  the  equation  which  the  stress  function  satisfies  are  the  three 
fourth-order  differential  terms  which  constitute  the  biharmonic  equation,  and  the  terms  in 
the  equation  involving  material  gradients  are  the  lower  order  differential  terms.  These  lead 
to  that  the  equations  for  >’)  and  p  are  the  same  as  those  for  homogeneous  materials. 
The  stress  intensity  factors  Kh  Ku  and  Km  are  functions  of  the  material  gradients,  external 
load,  and  geometry.  Material  gradients  do  not  affect  the  order  of  the  singularity  and  the 
angular  functions,  but  do  affect  the  stress  intensity  factors.  As  a  result  the  near-tip  stresses 
have  the  same  form  as  that  for  a  homogeneous  material.  For  an  interface  crack,  stresses 
have  an  oscillatory  singularity,  and  both  the  stress  intensity  factors  and  angular  functions 
involve  Dundurs’  parameters,  i.e.. 


M*=P(6]  +H-  62)+M 

Fig.  1 .  A  semi-infinite  crack  in  a  strip  of  a  functionally  graded  material  subjected  to  edge  loading. 
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Re(tfr*)  /n  x  ,  Im(AYc)  N  ,  Km  _ni,m 
<7,7  =  — 7=^  <7„  (0.  e)  +  — 7=-  <7,,  (0,  e)  +  — p=  a"1  (0), 
n/2tc/-  V-71'’  V-77'’ 

where  £  =  A^1  +  /A^11  is  complex  stress  intensity  factor,  and 


(8) 


1  1-jS 
e  =  —  In  - — - . 
2n  l+P 


(9) 


In  (9),  P  is  one  of  the  two  Dundurs’  parameters.  The  Dundurs'  parameters,  a  and  P,  are 
defined  as 


_  +  +  1) 

*  ^i(k'2  +  1)+/<2(ki  +  1)’ 

o  /<i(>v'2-1)-/<2(k'i -1) 

^,(K'2  +  l)  +  /i2(K,  +  l)’ 


(10) 


where  /q  and  /u2  are  the  shear  moduli  of  the  two  bulk  materials ;  k,  =  3  -  4v,  for  plane  strain 
and  K)  =  (3  —  v,)/l  +v,)  for  plane  stress  (/  =  1,  2),  with  v,  and  v2  being  the  Poisson’s  ratios 
of  the  two  bulk  materials.  It  is  noted  that,  by  considering  material  gradients  near  the  tip  of 
an  interface  crack,  the  oscillatory  behavior  is  removed,  and  the  angular  functions  become 
independent  of  material  properties.  In  this  sense,  the  solutions  presented  here  represent  a 
more  physically  acceptable  description  of  interface  crack  tip  fields,  at  least  for  interfaces 
that  have  a  finite  width. 

The  strains  obtained  from  the  stresses  given  in  (7)  are 

%  =  S,Jkl(0)aki  +  [Sijki(y)  -  Sijk,(0)]crkh  (11) 

where  SiJki(y )  is  the  compliance  tensor,  and  Sy*,( 0)  is  the  tensor  at  the  crack  tip.  The  second 
term  in  the  above  equation  is  in  the  order  of  rv:.  So  the  singular  strain  field  is 


£,7  =  S,jkl(0)akh  (12) 

From  (12),  one  is  able  to  show  that  the  near-tip  displacement  field  is  the  same  as  that  for 
the  homogeneous  materials. 

From  (7),  the  traction  at  the  distance  ;•  ahead  of  the  crack  tip  is 


K 

Oyy  +  K 7xr  = 

for  an  in-plane  problem.  For  a  mode  III  problem,  the  traction  is,  likewise, 


(13) 


Jlnr 


(14) 


Having  the  near  tip  stress  and  displacement  fields,  the  energy  release  rate  of  the  crack  tip 
is  obtained  as 


<$  = 


Kj 

£'(0) 


Kj ■  Kjn 

+  E\ 0)  +  2/2(0)  ’ 


(15) 


where  E'(0)  and  /<( 0)  are  the  Young’s  modulus  and  the  shear  modulus  at  the  crack  tip, 
respectively.  It  can  be  seen  that  the  above  eqns.  (7)  and  (12)-(15),  are  independent  of  the 
forms  of  the  material  properties  and  the  orientation  of  the  crack,  and  they  all  have  the 
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same  forms  as  those  for  homogeneous  materials.  The  path-independence  of  the  J  integral 
(Rice.  1968)  holds  if  the  crack  is  perpendicular  to  the  direction  along  which  material 
properties  change ;  this  is  implied  in  Rice’s  original  proof  for  homogeneous  materials.  Using 
the  near  tip  fields  obtained  above,  it  can  be  shown  that  the  J  integral  is  equal  to  the  energy 
release  rate  for  the  crack  perpendicular  to  the  direction  along  which  material  properties 
change. 

The  complex  stress  intensity  K  =  Ki  +  iKu  for  FGMs  has  the  following  generic  form 

K=  \K\e*9  (16) 

where 

\j/  =  tan-1  (17) 

Aj 


is  the  phase  angle  of  the  complex  stress  intensity  factor.  The  phase  angle  measures  mode 
mixity,  i.e.9  the  proportion  of  the  shear  traction  to  the  normal  traction  ahead  of  the  crack 
tip.  since 


\p  =  tan  1  (• — ^  .  (18) 

\a>y/e= o.r-o 

As  a  result  of  the  regular  singularity,  this,  again,  is  consistent  with  the  phase  angle  defined 
for  cracks  in  homogeneous  materials.  In  the  case  of  interface  cracks,  a  material  length  is 
needed  to  define  the  phase  angle. 

As  a  starting  point,  we  postulate  that  the  crack  starts  to  propagate  when  the  energy 
release  rate  reaches  a  critical  value  T,  the  toughness  of  the  FGM.  The  toughness  is  likely 
dependent  on  the  material  gradients,  the  position  of  the  crack  tip,  namely  h/H  for  the 
configuration  shown  in  Fig.  1,  and  the  mode  mixity  ij/.  It  is  also  possibly  dependent  on  the 
propagation  direction  </>(-  n  <  (j)  <  n)  which  is  the  angle  between  the  propagation  direction 
and  the  x  axis  in  Fig.  1.  The  energy  release  rate  is  a  function  of  the  external  load,  elastic 
constants,  the  angle  <p  and  the  mode  mixity  \j/.  Now,  the  fracture  criterion  is  stated  as 

^  =  r,  ^(sr-n  =  o.  O9) 


The  criterion  also  determines  the  propagation  direction  (kink  angle).  If  toughness  in  the 
direction  <j>  #  0  is  relatively  larger  than  that  in  the  direction  <f>  =  0,  the  crack  would 
propagate  along  its  original  orientation.  In  this  case,  the  fracture  criterion  is 


g  =  r 


(20) 


where  r(h/H,  i p)  is  the  toughness  along  the  direction  <f>  =  0.  The  toughness  of  FGMs 
may  be  measured  by  experiments  or  obtained  from  micromechanics  by  considering  their 
microstructures. 

4.  THE  IN-PLANE  PROBLEM 

The  in-plane  crack  problem  is  illustrated  in  Fig.  1,  where  a  semi-infinite  crack  in  a 
strip  occupies  the  negative  x  axis  and  the  crack  tip  is  at  the  origin.  The  material  properties 
change  along  the  y  axis.  The  geometry  is  specified  by  /?,  the  distance  between  the  crack  face 
and  the  upper  boundary,  and  H ,  the  distance  between  the  crack  face  and  the  lower  boundary. 
The  body  extends  infinitely  in  both  the  positive  and  negative  x  axes,  and  is  loaded  at  the 
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left  side  far  behind  the  crack  tip.  The  deformation  far  behind  the  crack  tip  consists  of 
bending  and  compression  in  the  upper  arm.  and  bending  and  tension  in  the  lower  arm.  The 
external  load  results  in  two  forces  which  are  of  the  same  magnitude  P  but  which  act  in 
opposite  directions,  and  two  bending  moments,  M  and  M*.  The  compressive  and  tensile 
forces  act  at  the  neutral  axes  of  the  upper  and  lower  arms,  respectively.  One  of  the 
two  moments,  M,  is  given  independently,  and  the  other  is  M*  =  +  77— <52)  by 

equilibrium,  where  <5i  is  the  distance  between  the  neutral  axis  of  the  upper  arm  and  the 
crack  face  and  S2  is  the  distance  between  the  neutral  axis  of  the  lower  arm  and  the  lower 
boundary.  According  to  the  small-strain  compatibility  equations,  the  only  non-zero  strain, 
exx,  in  the  two  arms  far  behind  the  crack  tip  varies  linearly  along  they  axis. 

We  consider  bending  deformation  in  the  two  arms  far  behind  the  crack  tip  resulting 
from  the  two  bending  moments.  The  strain  exx  can  be  expressed  as 

exx  =  -K'(>— 5,)  (21) 

in  the  upper  arm,  where  k  is  the  curvature  of  the  upper  arm.  The  moment  M  relates  to  the 
curvature  by 


M  =  k£0/,  .  (22) 

In  (22),  7,  is  the  moment  of  inertia  and  is  given  by 

r  h 3 

/,  =  I  e-'Xy-S,)2  d y  =  —a,  {yh),  (23) 

where 


«i(V*)  =  12 


yh 


2  <5,  2  ~ 

{yh?  h  (y/OT 

(24) 


When  y  —  0,  a,  =  1  so  that  7,  =  A3/ 12,  a  standard  result  for  homogeneous  materials.  The 
only  non-zero  stress  in  the  arm  is  the  normal  stress  in  the  cross  section, 


Miy-S? 
- ; - e,J 


(25) 


By  the  equilibrium  requirement  axxdy  =  0,  the  position  of  the  neutral  axis  is  obtained  as 


<3,  _  yh  e-7' — tyh  + 1 

h  ~  y/7(e'"- 1) 


Similarly,  results  for  the  lower  arm  are 


oxx  = - - -e-'\ 

h 


h  =  ^^h-m{r-S2?dy  =  ^^ 


12 


6  2  yH  e'w  -  + 1 

H  ~  yH(t/H  —  1 )  ’ 


(26) 


(27) 


where 
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Mr#)  =  . 
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(28) 

We  next  consider  axial  deformation  in  the  two  arms  resulting  from  the  two  forces.  Since 
the  compressive  force  acts  on  the  neutral  axis  of  the  upper  arm,  the  deformation  caused  by 
it  is  uniform  compression.  The  stress  distribution  in  the  cross  section  to  produce  the  uniform 
compression  is 


&xx  = 


(29) 


Similarly,  the  stress  distribution  in  the  cross  section  of  the  lower  arm  due  to  the  tensile 
force  is 


Py 

1  —  e“v// 


e7*v. 


(30) 


It  is  seen  that  in  the  cross  sections  the  strains  vary  linearly  with  y9  whereas  the  stresses  vary 
exponentially  with  y.  Expressions  for  the  stresses  and  strains  for  the  case  that  material 
properties  change  linearly  along  the  y  axis  have  been  derived  by  Freund  (1993)  and 
Giannakopoulos  et  ai  (1994)  in  studying  thermal  and  mechanical  responses  in  a  com- 
positionally  graded  layer  sandwiched  between  two  dissimilar  materials  without  cracks, 
where  a  linear  stress  distribution  along  the  thickness  of  the  compositionally  graded  layer  is 
obtained  from  the  linear  variation  of  material  properties.  Giannakopoulos  et  aL  (1994) 
have  also  investigated  plastic  deformation  in  the  compositionally  graded  layer.  More 
recently  Maewal  et  al  (1995)  have  developed  a  more  general  framework  for  analyzing 
thermally  induced  stresses  in  generally  orthotropic  FGMs  with  arbitrary  gradients. 

Having  the  remote  field,  the  complex  stress  intensity  factor  K  is  obtained  by  the 
application  of  the  J  integral,  dimension  analysis  and  linearity  consideration.  The  procedure 
is  similar  to  that  in  Suo  and  Hutchinson  (1990)  for  an  interface  crack.  The  complex  stress 
intensity  factor  is  obtained  as 


where 


K-Kl  +  tKa-j^fip-i'»J±M  )e- 


1  1 
A  =  yh[— — r  + 


-1  l-e-'V  y-2\H  h 


+-(4^+i-  SA!h 


H  H ’ 


,.12+lHf‘Y 

a,  «2  \HJ 

12  (h  <5,  ,  S2\(h\2 

in<P'«2yr/U  h)\h)  ■ 


(31) 


(32) 


In  (31),  a)  =  co(>7?,  h/H)  is  to  be  determined,  and  is  in  the  range  0  ^  co  ^  nil.  The  complex 
stress  intensity  factor  is  fully  obtained  apart  from  the  dimensionless  real  scalar  a>(yhJi/H). 
The  expression  for  K  has  a  similar  form  as  that  for  an  interface  crack.  In  the  interface  crack 
case,  Dundurs’  parameters  enter  the  solution  as  variables  of  A ,  /,  </>  and  co,  whereas  the 
material  constant  y  enters  the  solution  as  a  variable  of  those  parameters  in  our  case. 

To  determine  co(yh,h/H),  we  solve  the  full  boundary  value  problem  for  given  yh  and 
/?///,  using  the  integral  equation  method.  Integral  equations  for  this  problem  can  be  obtained 
by  Fourier  transforms.  The  numerical  procedure  is  to  distribute  dislocation  densities  simu¬ 
lating  opening  and  sliding  displacements  of  the  crack  in  terms  of  Chebyshev  polynomials, 
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Table  1 .  u)(yh.  h  H)  (in  degrees) 


V  h 


h'jH 

0.2 

0.4 

0.6 

0.8 

1.0 

1.2 

1.4 

1.6 

1.8 

2.0 

2.2 

2.4 

0.0 

53.8 

54.9 

56.0 

56.9 

57.5 

58.0 

58.3 

58.5 

58.6 

58.6 

58.5 

58.4 

0.1 

53.5 

54.8 

56.0 

56.9 

57.5 

58.0 

58.3 

58.5 

58.6 

58.6 

58.5 

58.4 

0.5 

51.6 

52.4 

53.2 

54.0 

54.7 

55.3 

55.9 

56.4 

56.8 

57.1 

57.4 

57.5 

1.0 

49.6 

50.1 

50.7 

51.2 

51.7 

52.2 

52.7 

53.1 

53.5 

53.8 

54.1 

54.4 

2.0 

47.0 

47.3 

47.6 

47.9 

48.2 

48.5 

48.8 

49.1 

49.4 

49.6 

49.9 

50.2 

10 

40.8 

40.9 

40.9 

41.0 

41.1 

41.2 

41.2 

41.3 

41.4 

41.5 

41.6 

41.7 

100 

39.0 

39.0 

39.0 

39.1 

39.1 

39.1 

39.1 

39.1 

39.1 

39.1 

39.1 

39.1 

Fig.  2.  Numerical  results  for  a )(yh,  hjH). 


Po 

J 

Po 

1 

-  -  H 

_h 

H 

7L 

f4 

Fig.  3.  A  four-point  bending  specimen. 


and  to  adjust  the  coefficients  of  these  polynomials  to  satisfy  the  integral  equations  which 
are  expressions  of  equilibrium  (see  Thouless  et  al 1987).  Table  1  gives  numerical  results 
for  co  when  0  ^  yh  <  2.4  and  hjH  —  0,  0.1,  0.5.  1,2,  10  and  100,  which  are  also  shown  in 
Fig.  2.  For  most  cases,  the  co  increases  as  yh  increases,  and  the  increase  is  larger  for  smaller 
hjH.  For  hjH  =  100,  the  numerical  solution  shows  little  change  of  co  when  yh  varies  between 
0  and  2.4. 

When  M  =  0,  the  phase  angle  \j/ =  co;  when  P  —  0  (double-cantilever  beam). 
\j/  =  co -f<p  —  90°.  A  four-point  bending  specimen  configuration  shown  in  Fig.  3  can  be 
reduced  to  the  present  problem  by  cutting  it  from  the  middle.  By  a  superposition  scheme 
(Suo  and  Hutchinson,  1990)  and  above  deformation  analysis,  the  force  P  and  the  moment 
M  are 
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P  J  Ch 

P  =  ~y  (y+H-S^e  '  dy, 

*  Jo 

M  =  ^y  (y + // -  <5 3 ) ( v — <5 ,  )  e”'  dy,  (33) 

*  Jo 

where 

rH+h 

/=  (>'-<53)2e''°-mdj.  (34) 

Jo 

In  (33)  and  (34),  S 3  is  the  neutral  axis  of  the  right  side  of  the  plate,  and  is  given  by  replacing 
h  by  H+h  in  (26).  Figure  4  shows  that  the  phase  angle  and  the  magnitude  of  the  complex 
stress  intensity  factor  for  the  four-point  bending  specimen.  The  phase  angle  varies  as  yh 
varies  between  0  and  2.4.  The  variation  is  larger  for  smaller  h/H,  and  is  quite  small  when 
h/H  increases  to  10.  The  magnitude  of  the  complex  stress  intensity  factor  increases  as  yh 
increases,  and  significantly  increases  as  h/H  increases.  Figure  5  shows  the  phase  angle  and 
the  magnitude  of  the  complex  stress  intensity  of  the  double-cantilever  beam.  In  this  figure, 


3* 

Fig.  4.  The  phase  angle  if/  and  the  magnitude  of  the  complex  stress  intensity  factor  \K\  for  four- 

point  bending  specimen. 
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(a) 


Fig.  5.  The  phase  angle  and  the  magnitude  of  the  complex  stress  intensity  factor  |AT|  for  double¬ 
cantilever  beam. 


the  phase  angle  increases  as  yh  or  h/H  increases.  Similar  to  Fig.  5(a),  the  variation  of  the 
phase  angle  for  yh  between  0  and  2.4  is  insignificant  when  h/H  increases  to  10. 

The  ratio  of  the  Young’s  moduli  of  the  two  material  phases  in  a  FGM,  £'  and  £,',  is 
usually  less  than  10.  From  (5), 


In  (E'u/E',) 
1  +H/h 


^  In  10  <2.4. 


(35) 


This  tells  us  that  the  numerical  results  for  yh  between  0  and  2.4  given  here  provide  a 
complete  solution  for  the  semi-infinite  crack  in  FGMs,  shown  in  Fig.  1. 

To  obtain  a  quantitative  feel  for  the  behavior  of  the  solution,  consider  the  double¬ 
cantilever  beam  (/.<?.,  P  =  0)  with  E'JE',  =  7.  If  h/H  =  1,  from  (35),  yh  *  1.  The  stress 
intensity  factors  are  Kth32/M  =  3.55  and  Kulr'  2/M  =  1.02,  whereas  for  a  homogeneous 
material  with  h/H  =  1,  they  are  3.46  and  0,  respectively.  If  h/H  =  0.1,  yh  %  0.18.  The  stress 
intensity  factors  are  K,h*-/M  =  1 .92  and  Kuh-  2/M  =  - 1 .34,  whereas  for  a  homogeneous 
material  with  h/H  =  0.1,  they  are  1.96  and  -1.47,  respectively.  If  h/H  =  10,  yh  *  1.8. 
The  stress  intensity  factors  are  Kthy2/M  =  63.23  and  Kuh:''2/M  =  50.71,  whereas  for  a 
homogeneous  material  with  h/H  =  10,  they  are  62.09  and  46.52,  respectively.  This  shows 
that  the  change  of  KH  is  larger  than  that  of  Kt  due  to  the  change  of  material  gradients.  Also, 
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the  stress  intensity  factors  at  E'JE)  =  7  are  larger  than  those  of  a  homogeneous  material 
except  the  case  for  hjH  =  0.1 . 

For  the  four-point  bending  configuration  with  E'u/E'i  =  7  and  h/H=  1,  the  stress 
intensity  factors  are  K)h3/2/(Pof)  =  2.21  and  Knlf  ‘/(Pol)  =  2.07;  they  are  1.73  and  1.50, 
respectively,  for  a  homogeneous  material.  If  hjH  =  0.1,  the  stress  intensity  factors  are 
K\hy2j(PJ)  =  0.05  and  Kuhx2/(P0l)  =  0.06 ;  they  are  0.03  and  0.03,  respectively,  for  a  homo¬ 
geneous  material.  If  /?///  =  10,  the  stress  intensity  factors  are  Kilr  I(PqI)  =  62.79  and 
KnhV2/(P0I)  =  51.21 ;  they  are  61.30  and  47.45,  respectively,  for  a  homogeneous  material. 
Similar  to  the  double-cantilever  beam,  the  change  of  Ku  is  larger  than  that  of  K}  due  to  the 
change  of  material  gradients ;  and  the  stress  intensity  factors  at  EuIEl  =  7  are  larger  than 
those  of  a  homogeneous  material. 


5.  SOLUTION  FOR  ORTHOTROPIC,  FUNCTIONALLY  GRADED  MATERIALS 

In  this  section,  we  consider  the  strip  in  Fig.  1  is  made  of  an  orthotropic,  functionally 
graded  material.  The  problem  is  solved  by  orthotropy  rescaling.  In  general,  there  are  three 
material  gradients  associated  with  two  Young’s  moduli  and  one  shear  modulus.  The  moduli 
can  be  written  in  following  forms 

E\  OO  —  E\Qt’1’ , 


E2(y)  =  E20  e7--\ 

Uni) ;)  = /•«(>')  er-,v>  (36) 

where  yu  y2  and  y3  are  material  constants;  £10  and  £20  are  E\ ( v)  and  E'2(y)  at  y  =  0, 
respectively;  /i , 2(y)  is  the  shear  modulus  in  the  x— y  plane.  For  plane  stress  problems,  £| 
(>•)  =  £,(>•)  and£'2(y)  =  £2(y),  with  £,(y)  and  £2(y)  being  Young’s  moduli  in  the  directions 
parallel  to  the  x  axis  and  the  y  axis,  respectively.  For  plane  strain  problems,  £', 
(>)  =  £,(j)/[l-vi3(y)v310>)]  and  £',(y)  =  £2(y)/[l -v23(y)v32(y)],  with  v13(>0  and  v31(>-), 
and  v23(j!)  and  v32(y)  being  four  Poisson’s  ratios  in  the  x-z  and  y-z  planes,  respectively. 
For  isotropic  materials,  the  elastic  properties  in  (36)  reduce  to  those  in  (1)  and  (2).  The 
variation  of  Poisson’s  ratios  can  also  be  written  in  the  exponential  form  similar  to  that  of 
the  Poisson’s  ratio  in  (1)  for  isotropic  materials. 

We  consider  a  special  set  of  elastic  properties,  which  is  given  by 


Vi  =  y  ~  =  y3  = 

Vi,  (A)  =  v2,o(l  +£.v)  evv, 


v  1  2  ( t  )  —  vl  20  (i  "b  ®3’)  ®”  » 


1 2  C.v) 


F2(v) 

2tVI+v;,(v)) 


(37) 


whereas  v210  and  v120  are  v'21(y)  and  v'12(y)  at  y  =  0,  respectively;  e  is  a  constant;  and 


£2o 


(38) 


For  plane  stress  problems,  v'21(y)  =  v2l(y)  and  v',2(y)  =  v,2(y),  with  v21(y)  and  v12(y) 
being  two  Poisson’s  ratios  in  the  x—y  plane.  For  plane  strain  problems, 
v',(y)  =  [v12(y)  +  v13(y)v32(y)][l -vI3(y)v31(y)],  and  v'21(y)  =  [v2I(y)  + v23(y)v31(y)]/ 
[1  —  v23(y)  v32(  v)] .  These  forms  for  the  material  properties  of  orthotropic,  functionally  graded 
materiais  provide  analytical  flexibility,  and  lead  to  somewhat  simple  forms  for  the  field 
equations. 
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Using  this  set  of  elastic  properties,  the  stress  function  <b(.\-,y)  satisfies  the  following 
equation 


C4Q  p;  d4®  d4® 


cx 


6x28y2 


8vA 


c"?®  .c5®  . .  d2<t> 

/.  2y/.  —  +  y2  /  —  =  0. 

cx-cy  cy  dy 


(39) 


Making  a  variable  change 


x 


'f, 


eqn  (39)  becomes 


o4®  d4® 

dc4  +2dFd? 


d2®  0  /c-Q  g2®\  d2® 


5/  ~2y^+^J 


+yv=0 


(40) 


(41) 


in  the  t-y  plane.  The  above  equation  shows  the  orthotropy  rescaling  (Suo  et  ai,  1991) 
works  for  nonhomogeneous  materials  which  obey  (36)  and  (37),  since  It  is  the  same  as  eqn 
(3)  for  isotropic  materials.  Stresses  in  terms  of  ®  (f ,  y)  are 


c*® 

;-“I/2<T.vv 

cy 


d2® 

d£d/ 


Stress  intensity  factors  are  expressed  as 


=  lim  , 

t-oo-ov  -^2’ 

7.~'8Wn  =  —  lim  . 

i-O.v-O  V  15 


(42) 


(43) 


According  to  above  analysis  and  the  solution  for  isotropic  materials  in  Section  4,  the 
stress  factors  to  the  orthotropic  problem  are  given  by 


(44) 


where  A,  I,  cp  and  co  are  the  same  as  those  in  the  solution  for  isotropic  materials  in  Section  4. 

From  above  expression,  K \  and  A),  are  those  of  the  isotropic  solution  modified  by 
multipliers  /3/s  and  )}  \  respectively.  The  phase  angle  of  the  orthotropic  problem,  ^orth, 
relates  to  that  of  the  corresponding  isotropic  problem,  i j/is,  by 


tan(</'orih)  =  z  1 4  tan(i /'is). 


(45) 


The  effects  of/,  on  the  crack-tip  field  depend  on  its  value.  When  /.  >  1,  the  stress  intensity 
factors  are  larger  than  those  of  the  isotropic  case,  and  the  phase  angle  is  smaller  than  that 
of  the  isotropic  case ;  and  when  /.  <  1,  the  stress  intensity  factors  are  smaller  than  those  of 
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Fig.  6.  The  effects  of  A  on  the  phase  angle  and  the  magnitude  of  the  complex  stress  intensity  factor 
|X|  for  four-point  bending  specimen. 


the  isotropic  case,  and  the  phase  angle  is  larger  than  that  of  the  isotropic  case.  For  the  four- 
point  bending  specimen  configuration,  these  effects  of  /.  are  shown  in  Fig.  6. 

6.  THE  ANTI-PLANE  PROBLEM 

The  cracked  strip  shown  in  Fig.  1  subjected  to  anti-plane  deformation  (mode  III 
problem)  is  considered  in  this  section.  For  the  mode  III  problem,  we  take  the  shear  modulus 
in  the  following  form 

/*(>’)  =  fio  e'1.  (46) 

where  /j0  is  the  shear  modulus  at  y  =  0  and  y  represents  the  material  gradient.  The  strip  is 
loaded  at  the  left  side  far  behind  the  crack  tip,  and  the  traction  is 

(7 V-.  =  —  o  \  e*  v  (47) 


for  x  —+  —  oo  and  y  >  0,  and 
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for  x  ■ 


<7ir  =  cr;  e  • 

—  oo  and  >•  <  0.  The  traction  produces  two  uniform  strains, 


(48) 


£  vr  =  - 


2^o 


for  x  —  oo  and  y  >  0,  and 


.  (49) 


£v-  = 


^0 


for  *  -*■  —  co  and  y  <  0.  From  the  equilibrium  requirement, 


Ch 


P  = 


<tXz  d y  = 


r  o 


cx-_  d y. 


(50) 


(51) 


and  (47)  and  (48),  we  have 


Py 


<7,  = 


C2  = 


ev/'-l 

Py 

1  —e~yi 


(52) 


Having  the  remote  field,  the  stress  intensity  factor  is  readily  obtained  from  the  J 
integral  as 


^n. 


(53) 


The  normalized  stress  intensity  factor  K1Uy/h/P ,  which  is  equal  to  the  second  square  root 
in  above  equation,  is  only  related  to  the  dimensionless  group,  yh  and  yH.  It  increases  as  y 
increases.  When  y  =  0,  the  stress  intensity  factor  recovers  the  solution  for  a  homogeneous 
material ;  when  y  =  oo  or  H  =  0,  it  is  unbounded.  A  plot  of  the  stress  intensity  fbctor  is 
shown  in  Fig.  7. 


Fig.  7.  The  stress  intensity  factor  Km  for  the  anti-plane  problem. 
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Fig.  8.  A  comparison  of  stress  intensity  factors  obtained  from  different  forms  of  the  shear  modulus. 


When  the  shear  modulus  varies  linearly  along  the  y  axis,  it  is  written  as 


,  .s  _  \L«-V-I  HuH+Hih 

^y)  l  (  tt  y  *  ' 


h  +  H‘ 


h  +  H 


(54) 


where  \xu  and  \it  are  shear  moduli  at  the  upper  and  lower  boundaries  of  the  strip,  respectively. 
The  stress  intensity  factor  for  the  linear  variation  of  the  shear  modulus  is 


P  I  2(r  +  rj)  T  rr}  +  2r  +  rj\ 
rrj  +  2r+rj\  n  r  +  2ij+l  )' 


(55) 


where  r  =  fiJni  and  rj  =  h/H.  A  comparison  of  the  stress  intensity  factor  obtained  from  the 
linear  variation  of  the  shear  modulus  with  that  from  the  exponential  variation  of  the  shear 
modulus  is  shown  in  Fig.  8  for  h/H  =  1.  The  results  show  that  the  difference  between  them 
is  quite  small,  less  than  5%  in  the  range  considered.  When  the  crack  moves  to  the  ceramic 
side,  the  difference  between  the  two  solutions  becomes  smaller,  and  when  it  moves  to  the 
metal  side  the  difference  becomes  larger.  The  difference  is  less  than  0.4%  for  h/H  ~  0.1 ; 
and  is  less  than  8%  for  h/H  =  10. 


7.  DISCUSSION 

A  complete  solution  to  a  semi-infinite  crack  in  a  strip  of  an  isotropic,  functionally 
graded  material  is  obtained.  It  is  shown  that  material  gradients  have  strong  effects  on  the 
stress  intensity  factors  and  the  phase  angle.  For  the  double-cantilever  beam,  the  mode  I 
stress  intensity  factor  is  3.55  and  the  mode  II  stress  intensity  factor  is  1.02,  when  the  crack 
is  at  the  middle  of  the  strip  (h/H  =  1)  and  the  ratio  of  the  Young’s  modulus  at  the  upper 
boundary  to  that  at  the  lower  boundary  is  7 ;  for  homogeneous  material  with  the  same 
geometry,  they  are  3.46  and  0,  respectively.  For  the  four-point  bending  specimen  con¬ 
figuration,  the  mode  I  stress  intensity  factor  is  2.21  and  the  mode  II  stress  intensity  factor 
is  2.07,  when  the  crack  is  at  the  middle  of  the  beam  (h/H  =  1)  and  the  ratio  of  the  Young’s 
modulus  at  the  upper  boundary  to  that  at  the  lower  boundary  is  7 ;  for  a  homogeneous 
material  with  the  same  geometry,  the  two  stress  intensity  factors  are  1.73  and  1.50,  respec¬ 
tively.  These  results  show  that  the  increase  of  the  mode  II  stress  intensity  factor  due  to  the 
increase  of  the  material  gradients  is  significant,  in  other  words,  the  mode  II  stress  intensity 
factor  plays  an  important  role  in  the  fracture  of  FGMs. 

The  solution  for  isotropic  materials  is  extended  to  orthotropic,  functionally  graded 
materials  by  orthotropy  rescaling.  The  effects  of  the  orthotropy  on  stress  intensity  factors 
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and  the  phase  angle  are  explicitly  shown  in  the  orthotropic  solution.  In  the  orthotropic 
solution,  since  we  assume  a  special  set  of  material  properties,  the  orthotropy  is  measured 
by  one  parameter,  the  ratio  of  the  Young's  modulus  in  the  direction  of  material  property 
variation  to  that  in  the  direction  perpendicular  to  the  above  property  variation  direction. 
For  general  orthotropic,  functionally  graded  materials,  there  are  other  parameters  in 
addition  to  the  ratio  for  characterizing  the  orthotropy.  However,  it  seems  that  the  ratio  is 
the  most  important  parameter  to  be  considered. 

The  crack  propagation  is  the  competition  between  the  driving  force,  the  energy  release 
rate,  and  the  toughness  of  the  material,  i.e..  a  crack  starts  to  extend  when  the  former  one 
exceeds  the  latter  one.  The  FGMs  are  expected  to  have  considerably  larger  toughness  than 
corresponding  bimaterials  because  there  are  no  large  weak  planes,  such  as  interfaces,  e.g ,  a 
layered  structure  with  compositionally  graded  interlayers  is  expected  to  have  a  lareer 
toughness  than  that  obtained  by  bonding  these  layers  with  sharp  interfaces.  On  the  other 
hand,  the  energy  release  rate  of  a  FGM  is  at  the  same  level  as  that  of  the  corresponding 
bimaterial.  Consider  a  bimaterial  which  has  the  same  configuration  as  the  FGM  shown  in 
Fig.  1 ,  above  the  x  axis  is  material  #1  with  Young’s  modulus  E'u  and  below  the  x  axis  is 
material  #2  with  Young’s  modulus  E).  Figure  9  shows  the  comparison  of  the  energy  release 
rate  of  the  FGM  with  that  of  the  bimaterial  for  the  double-cantilever  beam  when 
0  <  h/H  ^  2.  In  the  calculation,  E'u/E'i  =  7,  and  the  Poisson’s  ratios  of  the  two  bulk  materials 
forming  the  bimaterial  are  taken  to  be  0.3.  For  the  bimaterial,  the  two  Dundurs’  parameters 
are  a  =  0.75  and  /?  =  0.21.  Our  calculation  shows,  at  hlH  =  1,  the  energy  release  rate  is 
13.55  for  the  FGM,  whereas  it  is  16.51  for  the  bimaterial ;  at  h/H  =  0.1,  they  are  5.51  and 
4.61  for  the  FGM  and  the  bimaterial,  respectively;  at  h/H  =  10,  they  are  6559  and  6517 
for  the  FGM  and  the  bimaterial,  respectively.  When  the  crack  is  at  the  middle  of  the  plate 
(fi/H  =  1),  the  energy  release  rate  of  the  FGM  is  smaller  than  that  of  the  bimaterial ;  when 
the  crack  is  very  close  to  the  upper  or  lower  boundary,  the  former  one  is  larger  than  the 
latter  one.  But  in  any  case,  the  two  energy  release  rates  are  at  the  same  level.  This  fact 
reveals  one  of  the  advantages  of  using  FGMs,  i.e.,  FGMs  can  be  subjected  to  higher 
external  loads  than  corresponding  bimaterials. 

The  crack  propagation  direction  follows  different  criteria  for  different  kinds  of 
materials.  For  homogeneous  materials,  a  crack  propagates  along  the  direction  in  which  the 
mode  II  stress  intensity  factor  is  vanished,  and  the  toughness  is  independent  of  the  propa¬ 
gation  direction  and  the  mode  mixity.  For  bimaterials,  the  propagation  direction  of  an 
interface  crack  is  decided  by  the  driving  force  and  the  toughness  of  the  interface  and  the 
two  bulk  materials.  If  the  toughness  of  the  bulk  materials  is  relatively  large,  the  interface 
crack  would  extend  along  the  interface,  otherwise,  kinking  is  favored.  For  FGMs,  their 
toughness  is  likely  dependent  on  the  material  gradients,  the  position  of  the  crack  tip,  the 


Fig.  9.  A  comparison  of  the  energy  release  rate  of  the  FGM  with  that  of  the  corresponding 

bimaterial. 
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propagation  direction  and  the  mode  mixity.  From  a  continuum  point  of  view,  the  propa¬ 
gation  direction  is  the  direction  at  which  the  difference  of  the  energy  release  rate  and  the 
toughness  reaches  a  maximum  value,  as  discussed  at  the  end  of  Section  3.  For  complete 
understanding  of  the  fracture  behavior  of  FGMs  and  a  fully  rationalized  FGM  charac¬ 
terization,  experiments  need  to  be  carried  out  and  more  specimen  configurations  need  to 
be  calibrated. 

Acknowledgements — This  work  is  supported  by  the  Office  of  Naval  Research  through  Grant  N000 14-93-1-1 164. 


REFERENCES 

Delale.  F.  and  Erdogan.  F.  (1993).  The  crack  problem  for  a  nonhomogeneous  plane.  J.  Appl.  Mech.  50,  609-614. 

Delale.  F.  and  Erdogan.  F.  (1988).  On  the  mechanical  modeling  of  the  interfacial  region  in  bonded  half-planes. 
7.  Appl  Mech.  55,  317-324. 

Dundurs.  J.  (1969).  Edge-bonded  dissimilar  orthogonal  elastic  wedges.  J.  Appl .  Mech.  36,  650-652. 

Freund,  L.  B.  (1993).  The  stress  distribution  and  curvature  of  a  general  compositional!}'  graded  semiconductor 
layer.  7.  Crystal  Growth  132,  341-344. 

Giannakopoulos,  A.  E.,  Suresh,  S.,  Finot,  M.  and  Olsson,  M.  (1995).  Elastoplastic  analysis  of  thermal  cycling: 
layered  materials  with  compositional  gradients.  Acta  Metall.  43,  1335-1354. 

Holt,  J.  B.,  Koizumi,  M.,  Hirai,  T.  and  Munir,  Z.  A.  (eds)  (1993).  Functionally  gradient  materials.  Ceramic 
Transactions ,  Vol.  34,  The  American  Ceramic  Society,  Westerville,  OH. 

Jin,  Z.  and  Noda,  N.  (1994).  Crack-tip  singular  fields  in  nonhomogeneous  materials.  7.  Appl  Mech .  61,  738-740. 

Maewal,  A.,  Asaso,  R.  J.  and  Dao,  M.  (1995).  Residual  stresses  in  thin  film  structures  with  functionally  graded 
materials.  In  preparation. 

Noda,  N.  and  Jin,  Z.  (1993).  Thermal  stress  intensity  factors  for  a  crack  in  a  strip  of  a  functionally  gradient 
material,  Int.  7.  Solids  Structures  30,  1039-1056. 

Rice,  J.  R.  (1968).  A  path  independent  integral  and  approximate  analysis  of  strain  concentration  by  notches  and 
cracks.  7.  Appl  Mech.  35,  379-386. 

Suo,  Z.  and  Hutchinson,  J.  W.  (1990).  Interface  crack  between  two  elastic  layers.  Jnt.  7.  Fract.  43,  1-18. 

Suo,  Z.,  Bao,  G.,  Fan.  B.  and  Wang,  T.  C.  (1 991).  Orthotropy  rescaling  and  implications  for  fracture  in  composites. 
lnt .  7.  Solids  Structures  28,  235-248. 

Thouless,  M.  D.,  Evans,  A.  G.,  Ashby,  M.  F.  and  Hutchinson.  J.  W.  (1987).  The  edge  cracking  and  spalling  of 
brittle  plates.  Acta  Metall  35,  1333-1341. 

Yamanouchi,  M.,  Koizumi.  M.,  Hirai,  T.  and  Shiota,  1.  (eds)  (1990).  Proceedings  of  the  First  International 
Symposium  on  Functionally  Gradient  Materials,  Sendai,  Japan. 

Yang,  W.  and  Shih  C.  F.  (1994).  Fracture  along  an  interlayer,  lnt.  7.  Solids  Structures  31,  985-1002. 


> 


Pcrg\.mon 


PII:  S1359-6454(9-)00405-3 


Acta  mater.  Yoi.  45.  No.  8.  pp.  3265-3276.  1997 
£  1997  Acta  Metallurgies  Inc. 
Published  by  Elsevier  Science  Ltd.  All  rights  reserved 
Printed  in  Great  Britain 


1359-6454  97  S17.00  +  0.00 


A  MICROMECHANICAL  STUDY  OF  RESIDUAL  STRESSES 
IN  FUNCTIONALLY  GRADED  MATERIALS 

MING  DAO,  PEI  GU,  AKHILESH  MAEWAL  and  R.  J.  ASARO 

Department  of  Applied  Mechanics  and  Engineering  Sciences.  0411,  Universitv  of  California,  San  Dieeo, 

La  Jolla.  CA  92093,  U.S.A. 

(Received  11  June  1996;  accepted  27  November  1996) 

Abstract — A  physically  based  computational  micromechanics  model  is  developed  to  study  random  and 
discrete  microstructures  in  functionally  graded  materials  (FGMs).  The  influences  of  discrete 
microstructure  on  residual  stress  distributions  at  grain  size  level  are  examined  with  respect  to  material 
gradient  and  FGM  volume  percentage  (within  a  ceramic-FGM-metal  three-layer  structure).  Both 
thermoelastic  and  thermoplastic  deformation  are  considered,  and  the  plastic  behavior  of  metal  grains  is 
modeled  at  the  single  crystal  level  using  crystal  plasticity  theory.  The  results  are  compared  with  those 
obtained  using  a  continuous  model  which  does  not  consider  the  microstructural  randomness  and 
discreteness.  In  an  averaged  sense  both  the  micromechanics  model  and  the  continuous  model  give 
practically  the  same  macroscopic  stresses;  whereas  the  discrete  micromechanics  model  predicts  fairly  high 
residual  stress  concentrations  at  the  grain  size  level  (i.e.  higher  than  700  MPa  in  5-6  vol%  FGM  grains) 
with  only  a  3005C  temperature  drop  in  a  Ni-AhOj  FGM  system.  Statistical  analysis  shows  that  the 
residual  stress  concentrations  are  insensitive  to  material  gradient  and  FGM  volume  percentage.  The  need 
to  consider  microstructural  details  in  FGM  microstructures  is  evident.  The  results  obtained  provide  some 
insights  for  improving  the  reliability  of  FGMs  against  fracture  and  delamination.  ©  1997  Acta 
Metallurgica  Inc. 


1.  INTRODUCTION 

Functionally  graded  materials  (FGMs)  are  spatial 
composites  within  which  the  composition  of  each  of 
the  two  material  phases  that  form  the  FGMs  varies 
along  their  thickness  direction.  The  variation  is 
designed  to  be  tailorable  so  as  to  achieve  predeter¬ 
mined  responses  to  given  mechanical  and  thermal- 
mechanical  loadings.  Within  a  FGM,  the  different 
material  phases  have  different  functions.  In  a 
metal-ceramic  FGM,  the  metal-rich  side  is  placed  in 
the  region  where  mechanical  performance,  such  as 
toughness,  needs  to  be  stronger;  and  the  ceramic-rich 
side,  which  has  better  thermal  resistance,  is  exposed 
to  high  temperatures,  or  placed  in  the  region  where 
there  is  a  potentially  severe  temperature  variation. 
Also,  FGMs  can  reduce  the  thermal  mismatch  at  the 
interfaces  of  bimaterials  and,  therefore,  largely 
reduce  the  possibility  of  fracture  caused  by  the 
mismatch.  Applications  of  FGMs  include  aerospace, 
power  generation,  furnaces  and  others  where  strong 
material  performance,  especially  the  ability  to  resist 
thermal  shock,  is  required  or  expected. 

Material  gradients,  induced  by  the  change  in 
material  properties,  make  FGMs  different  in  behav¬ 
ior  from  homogeneous  materials  and  traditional 
composite  materials.  Over  the  past  few  years,  there 
have  been  a  number  of  works,  both  theoretical  and 
experimental,  to  study  the  responses  of  FGMs  to 
mechanical  and  thermal  loads  under  various  loading 
conditions,  for  various  geometries  and  in  various 


deformation  and  fracture  mechanisms,  including 
elastic  and  plastic  aspects  and  crack  propagation 
[1-1 1).  Most  of  the  previous  studies  above  focused  on 
the  continuous  approach  which  considers  that  the 
material  properties  change  continuously,  as  shown  in 
Fig.  1(b).  The  continuous  model  gives  correct 
solutions  to  such  problems  as  elastic  deformation  in 
the  ceramic-rich  side  and  plastic  deformation  in  the 
metal-rich  side,  when  the  scale  considered  is  much 
larger  than  that  of  the  grain  sizes  of  the  constituent 
material  phases.  It  also  gives  a  good  prediction  for 
damage  initiation  from  an  imperfection,  such  as  a 
void  or  crack,  when  the  size  of  the  imperfection  is 
much  larger  than  the  grain  size. 

The  microstructures  in  FGMs  are  discrete  and 
random  in  nature,  as  schematically  shown  in  Fig.  Ha). 
The  strongly  heterogeneous  microstructure  is  likely, 
at  least  possible,  to  cause  locally  concentrated 
residual  stresses  during  thermal  or  mechanical 
loading.  These  locally  concentrated  stresses,  es¬ 
pecially  those  high  in  tension,  may  act  to  initiate 
small  cracks  and  voids.  The  development  of  these 
small-scale  failures  may  lead  to  large-scale  failures 
and  result  in  the  fracture  of  the  whole  structure. 
Experiments  on  Si-C  FGM  by  Sohda  et  aL  [11] 
show-ed  that  the  porous  microstructure  has  a  much 
better  resistance  against  delamination  and  crack 
propagation  than  the  companion  dense  microstruc¬ 
ture,  where  the  latter  has  a  higher  level  of  local 
stresses.  Also  pointed  out  by  Finot  et  aL  [6],  to  study- 
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the  local  stress  distributions  and  concentrations 
within  the  FGM,  microstructural  details  such  as  the 
heterogeneous  microstructure  and  local  plastic 
deformation  must  be  considered  carefully.  The 
objective  of  this  study  is,  therefore,  to  explore  the 
microstructural  randomness  and  discreteness  vs 
macroscopic  material  gradients  and  geometries  with 
respect  to  the  thermal  residual  stresses  and  local 
residual  stress  concentrations  within  FGM  micro- 
structures.  A  discrete  computational  micromechanics 
model  is  developed.  In  our  discrete  micromechanics 
model,  the  ceramic  grains  are  treated  to  be  elastically 
deformed  as  the  typical  ceramic  materials;  the  metal 
grains  undergo  thermo-elastoplastic,  finite  defor¬ 
mation,  and  are  treated  using  crystal  plasticity 
theory.  The  results,  as  will  be  seen  in  the  later 
sections,  show  that  local  stress  concentration  at  the 
grain  size  level  is  significant.  For  the  purpose  of 
comparison,  we  also  solve  the  problem  by  the 
continuous  model. 

The  plan  of  the  paper  is  as  follows.  Both  the 
continuous  and  discrete  models  are  described  in 
Section  2,  where  a  brief  description  is  given  of  the 
crystal  plasticity  theory  used  for  metal  grains. 
Numerical  results  are  presented  in  Section  3.  In 
Section  3.1,  results  using  the  continuous  model  are 
presented,  where  influences  of  different  gradients  on 
macroscopic  residual  stresses  are  reviewed.  In  Section 
3.2.  results  using  the  discrete  micromechanics  model 
are  presented,  the  macroscopic  residual  stresses  as 
well  as  the  local  stress  concentrations  are  explored 
using  different  material  gradients  and  FGM  volume 

Material  Gradient 


percentages.  In  Section  3.3,  the  contribution  of 
plastic  deformation  within  the  discrete  micrcmechan- 
ics  mode!  is  studied.  In  Section  3.4.  a  short  summary 
is  given  on  the  statistical  analysis  of  the  residual  stress 
concentrations  with  respect  to  material  gradient. 
FGM  volume  percentage  as  well  as  the  plastic 
relaxation.  Finally,  discussions  and  conclusions 
follow  in  Section  4. 


2.  THE  CONTINUOUS  AND  DISCRETE  MODELS 

The  model  geometry,  as  shown  in  Fig.  2,  consists 
of  three  layers:  the  ceramic  layer  is  on  the  left  side; 
the  metal  layer  is  on  the  right  side:  and  the  FGM  is 
sandwiched  between  them.  As  a  model  system,  we 
choose  the  metal  to  be  Ni,  and  the  ceramic  to  be 
A1:03  in  this  study.  The  FGM  is,  therefore,  Ni-Al203 
FGM.  Both  continuous  and  discrete  models  for  the 
FGM,  including  the  numerical  consideration,  are 
described  below. 

2.1.  The  continuous  model 

We  define  .v  as  the  relative  distance  from  the 
ceramic-FGM  interface,  i.e.  a*  =  0  stands  for  the 
ceramic-FGM  interface  and  .v  =  1  stands  for  the 
FGM-metal  interface.  For  the  continuous  model  for 
the  FGM,  the  effective  material  properties  are 
assumed  to  follow  the  “rule  of  mixture'’: 

/4(a)  \  Melat(-V)/4 Mcul  4*  Fceramic(-V)^ Ceramic,  (1) 

Material  Gradient 


□  Ceramic  Grain 
01  Metal  Grain 


Material 

Property 


(a)  Discrete  Microstructure  (b)  Continuous  Model 


Fig.  1.  Schematic  drawings  of  functionally  graded  materials  (FGMs):  (a)  discrete  and  random 
microstructure  in  reality,  and  (b)  continuous  gradient  modeling  often  used. 


where  A  stands  for  either  the  elastic  constants,  E 
(Youngs  modulus)  and  v  (Poisson's  ratio),  or  the 
thermal  expansion  coefficient  «;  I'm^iCy)  and 
Ceramic  (a*)  are  the  volume  fractions  of  metal  and 
ceramic,  respectively,  at  the  position  .v.  The  simplified 
materia!  property  form  overlooks  the  interactions  of 
the  two  material  phases  at  the  microscopic  level,  so 
it  leads  to  an  approximate  solution.  The  more 
accurate  material  property  variation  form  at  the 
macroscopic  level  requires  a  better  understanding  of 
FGM  microstructure  and  its  deformation,  which  are 
the  focus  of  this  study.  We  will  only  obtain  the  elastic 
solution  for  the  continuous  model,  and  it  is  mainly 
for  comparison  with  the  solution  obtained  by  the 


discrete  model.  Plastic  deformation  of  the  sandwich 
structure  was  studied  using  a  continuous  model  in 
Giannakopoulos  et  al.  [5]  and  Finot  et  al .  [6]. 

The  thermoelastic  solution  in  this  case  may  be 
obtained  analytically  [5].  Here,  a  finite  element 
method  is  used  for  the  purpose  of  examining  its 
accuracy  for  FGMs.  In  the  implementation,  the 
FGM  layer  is  divided  into  30  micro-layers,  as  shown 
schematically  in  Fig.  2(a).  and  the  material  properties 
of  each  micro-layer  are  taken  to  be  constants. 

2.2.  The  discrete  mkromechanics  model 

We  have  developed  a  computational  micromechan¬ 
ics  model  for  FGMs  using  the  crystal  plasticity 
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theory.  Figure  2(b)  shows  the  model  geometry  for  a 
FGM  consisting  of  the  metallic  and  ceramic  grains 
randomly  distributed  within  it.  The  macroscopic 
material  properties  obtained  by  statistical  processing 
of  the  random  distribution  of  the  metal  and  ceramic 
grains  vary  continuously  along  the  thickness  direc¬ 
tion.  and  give  their  desired  variation  forms.  Each  of 
the  metal  grains  has  its  own  crystal  orientation  (also 
randomly  distributed)  which  is  shown  by  the  angle  \j/ 
in  Fig.  3,  and  its  thermoplastic  behavior  is  assumed 
to  be  governed  by  crystal  plasticity  theory.  The 
ceramic  grains  are  modeled  using  the  standard  linear 
elasticity  theory. 

The  two-dimensional  idealization  shown  in  Fig.  3 
was  introduced  by  Harren  el  al.  [12).  Harren  and 
Asaro  [13]  and  McHugh  ei  al.  [14]  for  f.c.c.  or  b.c.c. 
polycrvstals  and  their  metal  matrix  composites.  The 
three  slip  systems  are  arranged  in  an  equilateral 
triangle,  and  the  reference  laboratory  base  vectors  X, 
are  at  an  angle  ip  with  respect  to  reference  crystal  base 
vectors  a;.  The  slip  directions  in  this  model  geometry, 
Si,  s:  and  s?,  represent  the  close-packed  directions  of 
an  assemblage  of  close-packed  circular  cylinders. 
Since,  in  a  two-dimensional  model  two  independent 
slip  systems  can  accommodate  arbitrary  increment  of 
plastic  strain,  the  three  independent  slip  systems  here 
resemble  the  redundancy  exhibited  by  both  f.c.c.  and 
b.c.c.  crystals.  We  note  that  using  traditional  metal 
plasticity  theories  (i.e.  J2  flow  theory)  would  give  us 
similar  results  for  Ni  (f.c.c.).  If  any  low-symmetry 
crystals  (say  NiAl  or  TiAl)  are  involved,  then  crystal 
plasticity  theory  is  necessary  to  account  for  the 
orientation  dependent  deformation  behavior. 

The  single  crystal  constitutive  theory,  in  its  present 
form,  was  developed  by  Asaro  and  his  coworkers 
[14-21].  The  theory  which  will  be  briefly  described 
below  builds  on  the  pioneering  work  by  Taylor  [22] 
and  Hill  and  Rice  [23]. 


X, 


Fig.  3.  Two-dimensional  model  single  crystal  slip  geometry 
used  for  metal  grains.  The  three  slip  systems  are  arranged 
in  an  equilateral  triangle,  and  the  reference  laboratory  base 
vectors  X  are  at  an  angle  i, ft  with  respect  to  reference  crystal 
base  vectors  a,. 


The  total  deformation  gradient  is  decomposed  into 
plastic  (Fp).  thermal  (F'),  and  lattice  (F*)  pans,  as 
shown  in  the  insert  of  Fig.  2(b).  If  u  is  the 
displacement  vector  and  X  the  material  position 
vector  with  respect  to  the  reference  (undeformed) 
state,  F  =  I  +  cu  cX  (I  is  the  second-order  identity 
tensor)  and 

F  =  F*-F”-Fp.  (2) 

Plastic  deformation  occurs  by  the  flow  of  material 
through  the  lattice,  via  simple  shearing,  across  planes 
with  unit  normals  m,  and  in  directions  s*:  here  m7  and 
s,  represent  a  crystallographic  slip  plane  normal  and 
a  slip  direction,  respectively,  and  a  is  an  index  that 
designates  a  slip  system.  If  y,  is  the  slip  rate  on  the 
y.  slip  system,  then  the  velocity  gradients  of  this 
plastic  shear  flow  can  be  written  as 

FpFp-'  =Iv,s,m„  (3) 


where  the  summation  is  over  all  active  slip  systems. 
The  thermal  parts  of  the  velocity  gradients  are 
described  as 


F!j-F,_  1  =  $a;  a  =  J  ^  aya,a/,  (4) 

i  j 

where  6  represents  temperature  and  a  is  a  tensor 
whose  components,  ay,  with  respect  to  the  time 
independent  Cartesian  base  vectors,  a,,  are  the 
thermal  expansion  coefficients.  The  base  vectors  are 
aligned  with  the  crystal  lattice  in  the  reference 
configuration  in  some  standard  way,  e.g.  in  cubic 
crystals.  It  is  most  convenient  to  align  the  a;  base 
vectors  with  the  cube  axes,  in  which  case  a  would  be 
diagonal  with  all  components  equal. 

In  general,  y,  will  be  a  function  of  temperature, 
stress  state  and  material  state.  As  a  specific  example 
we  have  used  expressions  such  as 


‘A  =  a  sgn{ra)  j 


(5) 


to  represent  strain  rate  sensitive  power-law  type 
behavior,  In  the  expression,  n  is  the  material  rate 
sensitivity  parameter,  t*  is  the  resolved  shear  stress  on 
the  slip  system  a,  and  g,  >  0  is  its  current  strength. 

The  slip  system  hardness,  g,,  is  obtained  by  the 
path-dependent  integration  of  the  evolution 
equations  of  the  form: 


gi  —  X  ^L./?(7j)IV/il  +  gt6: 


IlV.ld/.  (6) 
*0  * 


where  ya  is  the  accumulated  sum  of  slips,  h,g  is  a 
matrix  of  hardening  moduli  and  gl  is  the  rate  of 
change  of  slip  system  hardness  with  respect  to 
temperature  alone.  The  initial  condition  for  this 
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Table  ).  Elastic  propenies  used  for  the  computations 


£  (GPa) 

r 

7  (K-‘) 

A1:0; 

3S0 

0.25 

7.4  x  10-* 

Ni 

214 

0.35 

15.4  x  I0-* 

evolution  is  given  by 

*,('/.  =  0, 6)  =  go(0),  (7) 

where  0  is  the  temperature. 

The  detailed  development  of  these  constitutive 
relations  can  be  found  in  McHugh  ex  al.  [14).  This 
constitutive  theory  has  been  implemented  into  finite 
element  codes,  using  a  rate  tangent  method 
introduced  by  Pierce  el  al.  [18]. 

3.  RESULTS 

The  geometry  of  the  sandwich  structure  is  specified 
in  Fig.  2.  To  avoid  the  edge  effects,  the  total 
length*/ width  ratio  (U/{H  +  2/;))  was  chosen  to  be  5, 
while  only  the  center  part  with  length  L(L  <  0.25Lo) 
was  considered  when  performing  residual  stress 
analysis.  Two  sets  of  FGM  interlayer  thickness,  i.e. 
H/(H  +  2 h)  =  40  and  70  vol%,  were  used. 

To  make  this  complicated  boundary  value  problem 
manageable  computationally,  a  relatively  coarse 
mesh  (i.e.  four  triangular  elements  per  grain)  is  used 
in  this  study  and  simple  square  grains  are  employed. 
Doing  so,  as  Taylor  [22,  24]  and  many  others  [25-27] 
did  successfully  in  modeling  polycrystalline  materials, 
effectively  treats  the  deformation  within  each 
individual  grain  as  uniform.  The  finite  element  model, 
however,  takes  the  interactions  between  all  constitu¬ 
ent  grains  into  consideration,  which  is  not  achievable 
using  Taylor  or  Sachs  type  models.  In  studying  the 
interactions  at  the  grain  size  level,  this  model  design 
is  at  least  a  first-order  approximation.  Also,  in 
keeping  all  the  microstructural  “building  blocks”  (i.e. 
ceramic  and  metal  grains  in  the  discrete  model) 
exactly  the  same,  the  relative  importance  of  material 
gradient  and  FGM  volume  percentage  can  be 
identified. 

The  specific  FGM  system  is  the  Al:Oj-Ni  system. 
In  this  study,  the  major  focus  is  on  the  relative 
importance  of  the  discrete  and  random  microstruc¬ 
ture  vs  FGM  volume  fraction  and  gradient  functional 
form.  With  that  in  mind,  and  noting  that  all  the  case 
studies  shown  later  were  performed  under  a  relatively 
small  (300:C)  temperature  variation,  for  simplicity 
the  elastic  properties  were  taken  to  be  constants.  The 
material  properties  used  for  the  computations  are 
listed  in  Tables  1  and  2.  A  two-dimensional  plane 


Fig.  4.  Different  gradient  functions  used  in  the  residual 
stress  analysis. 


strain  deformation  condition  was  imposed.  The 
thermal  loading  was  induced  by  cooling  the  structure 
by  30CFC.  The  temperature  was  assumed  to  be 
uniform  within  the  three-layer  structure,  and  the 
sandwich  structure  was  set  to  be  stress  free  at  the 
beginning  temperature.  For  the  metal  grains,  linear 
interpolation  was  used  to  obtain  the  temperature-de¬ 
pendent  critical  resolved  shear  stresses  go(0)  in 
equation  (7)  between  the  several  temperatures  shown 
in  Table  2.  The  shear  strain  hardening  h  in  equation 
(6)  was  taken  to  be  77.4  MPa,  which  was  convened 
from  the  data  in  Ref.  [28]  using  a  Taylor  factor  of 
3.06  for  f.c.c.  polycrystals  with  the  linear  hardening 
assumption.  A  low  material  rate  sensitivity  parameter 
is  given  as  n  =  0.005  in  equation  (5). 

We  write  KFgm  as  the  volume  fraction  of  Ni  within 
the  FGM  layer,  and  a*  as  the  relative  distance  from 
the  A1:0:-FGM  interface  (i.e.  a*  =  0  stands  for  the 
A1:03-FGM  interface  and  a*  =  1  stands  for  the 
FGM-Ni  interface).  As  shown  in  Fig.  4,  three 
functional  forms  were  used  in  the  computations: 

Linear:  I  fcm(.y)  =  (8a) 

FGM/1:  Ffcm(a')  =  1  -  (1  -  .v)"\  m  =  2,  4,  (8b) 

FGM/2:  KFOm(.v)  =  .v",  m  =  2.  4.  (Sc) 

When  m  =  1,  both  FGM/1  and  FGM/2  reduce  to  the 
linear  case  of  equation  (8a).  and  when  m  ^  2 
functions  in  the  FGM/1  class  have  zero  slope  at  the 
FGM-Ni  interface  (.v  =  1).  whereas  the  functions  in 
the  class  FGM/2  have  zero  slope  at  the  ALCR-FGM 
interface  (a*  =  0). 


T{  C) 

<r/  (MPa) 
tcr$s*  (MPa) 


Table  2.  Plastic  propenies  of  Ni  grains  used  in  the  computations _ _ 

:0  127  227  327  427  527  627  727  827 

148  153  140  13S  115  100  69  59  45 

48  4  50.0  45.8  45.1  37.6  32.7  22.5  19.3  14.7 


■'Data  from  Suresh  <■/  al.  [28]. 

'’Calculated  from  o.  with  a  Taylor  factor  of  3.06. 
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3.}.  The  continuous  model — macroscopic  stresses  is 
material  gradients  and  FGM  volume  percentages 

For  this  problem,  the  only  non-zero  stress 
component  is  the  in-plane  normal  stress  along  the  X: 
direction.  Figures  5(a)  and  (b)  show  the  results  of  the 
in-plane  normal  stress  distribution  along  the  A', 
direction  for  the  40  and  70  vol%  FGM.  respectively. 
Due  to  the  presence  of  the  FGM  layer,  the  stresses 
are  continuous,  and  have  continuous  derivatives 
wherever  the  material  property  variation  function  in 
equations  (8)  is.  It  is  noteworthy  that  although  the 
existence  of  the  FGM  layer  in  general  decreases  the 
stress  at  one  or  both  of  the  interfaces,  for  power  index 
m  ^  2  in  equations  (8),  there  is  an  extremum  in  the 
FGM  layer,  and  this  extremum  in  some  cases  has  a 
magnitude  close  to  that  of  the  normal  stress  at  the 
interface  in  the  base-line  case,  i.e.  a  sharp 
ceramic-metal  interface  without  the  FGM  layer  [29]. 
Comparing  Fig.  5(b)  with  Fig.  5(a),  it  is  found  that, 
in  general,  increasing  the  relative  FGM  volume 
percentage  decreases  the  stress  at  one  or  both  of  the 
interfaces.  These  results  are  similar  to  those  found  in 
Giannakopoulos  et  al.  [5]. 

Next,  we  will  explore  how  the  local  stress 
concentrations  interact  with  material  gradient  as  well 
as  the  FGM  volume  percentage.  The  averaged 
physical  meaning  of  the  continuous  solution  will  be 
clearer  when  we  present  the  discrete  solution  below 
and  compare  the  two  solutions. 

3.2.  The  discrete  model — local  stress  concentrations 
and  macroscopic  stresses 

In  this  section,  the  discrete  micromechanics  model 
is  used  and  only  elastic  deformation  is  considered. 
Plasticity  effects  will  be  studied  in  the  next  section. 

Figure  6  shows  contour  plots  of  (a)  <7::,  and  (b) 
averaged  in-plane  principal  stress  (p  =  (ou  +  0:2)1 2) 
developed  in  the  40  vol%  FGM  with  linear  gradient. 
It  is  clearly  seen  that  the  local  stress  concentration  is 
quite  high  and  the  stress  field  is  very  inhomogeneous, 
i.e.  the  stress  variations  among  many  of  the  adjacent 
grains  are  significant.  The  Gu  was  also  found  to  be 


inhomogeneous.  Due  to  the  thermal  mismatch 
between  ceramic  grains  and  metal  grains,  most  metal 
grains  experience  large  tensile  stresses;  and  this  is 
especially  true  for  metal  grains  near  the  ceramic- 
FGM  interface  and  those  in  the  middle  region  of  the 
FGM  layer. 

Additional  compulations  were  performed  using 
different  gradient  functions.  Figure  7  shows  contour 
plots  of  averaged  in-lane  principal  stress 
( p  =  (<7m  +  G22)/ 2)  developed  in  the  40  vol%  FGM 
with  (a)  gradient  function  EVcm  =  1  -  (1  -  .y): 
(FGM /l  m  =  2),  and  (b)  gradient  function  KFgm  —  .v2 
(FGM/2  m  =  2).  The  distribution  of  stresses  is  quite 
different  with  different  gradient  functions,  as  can  be 
seen  clearly  from  Figs  6(b).  7(a)  and  7(b).  Detailed 
examination  of  Figs  6  and  7  tells  us  that,  in  almost 
every  local  region  (say  take  5x5  grains  as  the  region 
size)  where  the  ceramic  grains  are  more  than  40  vol%, 
there  are  always  some  metal  grains  experiencing 
significant  tensile  stresses  for  all  three  cases.  The 
results  for  the  70  vol%  FGM  were  similar  to  those  for 
the  40vol%  FGM. 

For  such  discrete  . and  random  microstructures, 
more  physical  insights  can  be  gained  via  the  statistical 
analyses  of  the  stress  distribution.  Figure  8(a)  shows 
the  distribution  profiles  of  p  —  (an  4-  <7;;)/2  devel¬ 
oped  in  the  FGM  layer  (40  vol%  FGM)  with  three 
different  gradient  functions;  Fig.  8(b)  shows  the 
distribution  profiles  of  p  developed  in  the  FGM  layer 
(linear  gradient)  of  a  40  and  a  70  vol%  FGM, 
respectively.  The  three  distribution  profiles  in  Fig. 
8(a)  are  distinctively  different:  (i)  the  distribution 
peak  between  50  and  200  MPa  (mostly  in  metal 
grains)  drops  as  the  total  metal  composition  in  the 
FGM  layer  decreases;  and  (ii)  the  distribution  peak 
between  -50  and  -250  MPa  (mostly  developed  in 
ceramic  grains)  increases  as  the  total  ceramic 
composition  in  the  FGM  layer  increases.  The  FGM 
volume  percentage  also  has  large  influences  on  the 
distribution  profile  of  /;,  see  Fig.  8(b):  the  distribution 
peak  for  tensile  stresses  shifted  from  around  110  MPa 
(40  vol%  FGM)  to  about  50  MPa  (70  vol%  FGM); 
and  the  distribution  peak  for  compressive  stresses 


Fig.  5.  ln-p!ane  normal  stress  distributions  along  the  A'i  direction  using  the  continuous  model  with  (a) 

40  and  (b)  70  voi%  FGM,  respectively. 
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Fig.  6.  Contour  plots  of  (a)  (72:,  and  (b)  averaged  in-plane  principal  stress  (p  «  (ffn  +  <722 )/2)  developed 
within  the  linear  gradient,  40  vol%  FGM.  Only  elastic  deformation  is  considered  here. 


shifted  from  around  —150  MPa  (40vol%  FGM)  to 
about  —50  MPa  (70  vol%  FGM).  However,  the 
distribution  profile  for  high  stresses,  i.e.  for 
p  <  500  MPa  or  p  ^  —  500  MPa,  are  found  to  be_ 
insensitive  to  material  gradient  and  FGM  volume 
percentage. 

Finally,  for  the  purpose  of  comparison,  we 
averaged  the  stresses  over  each  column  of  elements  to 
get  the  mean  stress  along  the  vertical  direction. 
Figure  9  shows  the  macroscopically  averaged 
in-plane  normal  stress  developed  within  the  discrete 
microstructures  (shown  by  separated  symbols)  as 
compared  to  the  curves  obtained  using  the  continu¬ 
ous  model  (shown  by  continuous  lines).  It  is 
interesting  to  see  that,  although  there  are  a  lot  of  local 
stress  concentrations,  the  two  types  of  modeling 


approach  given  practically  the  same  averaged  (or 
macroscopic)  stresses.  In  an  averaged  sense,  only 
small  variations  can  be  found  in  Fig.  9  for  the  discrete 
microstructure  model  due  to  the  local  randomness 
and  discreteness.  This  shows  that  the  macroscopic 
stresses  are  those  based  on  a  scale  much  larger  than 
the  grain  size  and  obtained  without  considering 
certain  details,  such  as  local  stress  concentrations,  at 
a  smaller  scale. 

3.3.  The  discrete  model — influence  of  plastic  defor¬ 
mation 

In  this  section,  we  explore  the  effects  of  plastic 
deformation  within  the  metal  grains  in  the  discrete 
microstructure.  Besides  the  concentrated  stresses, 
large,  locally  concentrated  plastic  strain  accum 
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lation  during  repeated  thermal  cycling  may  also 
initiate  failure. 

Figure  10  shows  the  contour  plots  of  (a)  averaged 
in-plane  principal  stress  (p  =  (on  +  <?::)/ 2),  and  (b) 
accumulated  sum  of  slips  (y,)  developed  in  the  linear 
gradient,  40  vol%  FGM.  The  temperature  drop  was 
from  700  to  400°C,  and  the  plasticity  parameters  were 
taken  to  be  as  listed  in  Table  2.  Similar  to  the 
thermoelastic  case,  the  stress  distribution  is  again 
inhomogeneous,  with  many  metal  grains  experiencing 
high  tensile  stresses  and  many  ceramic  grains 
experiencing  high  compressive  stresses.  If  we  compare 
Fig.  10(a)  with  Fig.  6(b)  (elasticity  only),  the  stress 
concentration  in  the  ceramic  grains  is  significantly 
reduced  in  the  region  where  metal  content  is  greater 
than  70  vol%,  due  to  plastic  relaxation.  With  only  a 


300°C  temperature  drop,  Fig.  10(b)  shows  that:  (i) 
there  are  plastic  strain  accumulations  in  many  of  the 
metal  grains,  and  (ii)  certain  sites  have  relatively  high 
strain  accumulations,  about  1.5%.  The  high  strain 
accumulation  sites  seem  to  appear  in  the  regions 
where  metal  content  is  between  50  and  75  vol%. 

Figure  11  shows  the  distribution  profiles  of  p 
(  =  (ff,i  +  ff«)/2)  developed  within  the  FGM  layer 
(linear  gradient,  40  vol%  FGM),  for  both  the  elastic 
and  the  elastoplastic  case.  The  plastic  relaxation 
effect  is  very  clear  here  in  this  case,  where  stresses  are 
in  general  shifting  to  lower  magnitudes  with 
plasticity.  The  distribution  profile  for  high  tensile 
stresses  with  p  >  700  MPa,  however,  has  reduced 
only  slightly  with  plastic  relaxation.  Similar  to  the 
thermoelastic  case,  the  stress  distribution  profile  for 
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Fig  7  Contour  plots  of  averaged  in-plane  principal  stress  (p  =  (<xn  +  ou)fl)  developed  in  the  40  vol% 
FGM  with  (a)  gradient  function  Kfgm  =  1  -  (1  -  x)2  (FGM/1  m  «  2),  and  (b)  gradient  function 
Vfgm  =  x2  (FGM/2  m  =  2).  Only  elastic  deformation  is  considered  here. 
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Fig.  8.  (a)  Averaged  in-plane  principal  stress  (p  =  (<7U  +  on)l2)  distribution  profiles  developed  within  the 
FGM  layer  (40vo!%  FGM)  with  three  different  gradient  functions;  and  (b)  distribution  profiles  of  p 
developed  in  the  FGM  layer  (linear  gradient)  of  a  40  and  a  70  vol%  FGM,  respectively. 


high  tensile  stress  regions  is  insensitive  to  material 
gradient  and  FGM  volume  percentage.  On  the  other 
hand,  the  distribution  profile  for  high  compressive 
stresses  with  p  ^  —  250  MPa  (mostly  in  ceramic 
grains)  drops  significantly.  This  suggest  that,  when 
ceramic  grains  are  subject  to  tensile  stresses  if 
temperature  increases ,  the  plastic  relaxation  effects 
may  reduce  their  tensile  stress  concentrations. 

We  average  the  stresses  over  each  column  of 
elements  for  the  plastic  solution  to  obtain  the  averaged 
in-plane  normal  stress.  Figure  12  plots  the  macro- 
scopically  averaged  stress  for  both  the  elastic  and  the 
elastoplastic  case.  Comparing  Fig.  12  with  Fig.  10(b), 
it  is  found  that  the  metal  rich  section  and  part  of  the 
pure  metal  region  are  under  general  macroscopic 
yielding,  which  sets  the  maximum  magnitude  of  the 
macroscopic  stresses  for  the  plastic  case. 

3.4.  The  statistical  analysis  of  residual  stress  concen¬ 
trations — averaged  peak  stress 

Since  the  stress  distribution  profiles  of  p  for  the 
high  stress  area  shown  in  Figs  8  and  1 1  are  small,  to 


Relative  Distance 


Fig.  9.  Macroscopically  averaged  in-plane  normal  stress 
developed  within  the  discrete  microstructures  (shown  by 
separated  symbols)  as  compared  to  the  curves  obtained 
using  the  continuous  model  (shown  by  continuous  lines). 


get  more  reliable  statistical  results  we  employ  6  vol% 
APSP  to  treat  the  data  in  these  figures.  The  6vol% 
APSP  (averaged  peak  stress  of  p)  is  the  stress  p 
averaged  over  the  6  vol%  microstructure  of  the  FGM 
layer  which  has  the  highest  tensile  stresses  (p ). 
Similarly,  we  can  obtain  3  vol%  APSP,  9  vol%  APSP 
etc.  Figure  13  shows  the  6vol%  APSP  (averaged 
peak  stress  of  p)  for  different  material  gradients  and 
different  FGM  volume  percentages,  and  for  both  the 
elastic  solution  (marked  with  El)  and  the  elastoplastic 
solution  (marked  with  PI).  From  Fig.  13,  the 
distribution  profile  for  high  tensile  stresses  is  again 
found  to  be  insensitive  to  material  gradient  and  FGM 
volume  percentage,  and  the  plastic  relaxation  effect  is 
relatively  small  for  the  high  stress  regime.  Similar 
conclusions  can  be  drawn  from  the  3  vol%  APSP 
(averaged  peak  stress  of  p). 

The  above  statistical  results  are  interpreted  by  the 
following  observations.  As  mentioned  before,  high 
tensile  stresses  always  occur  in  those  metal  grains 
surrounded  by  many  ceramic  grains.  Comparing 
those  stresses  in  Figs  5  and  9,  the  stress  concen¬ 
trations  in  these  metal  grains  are  found  to  be  an  order 
of  magnitude  higher  than  the  macroscopic  stresses  in 
the  FGM  layer.  Considering  the  large  number  of  stiff 
ceramics  around  each  of  those  metal  grains,  the  high 
local  stresses  are  believed  to  be  insensitive  to  the 
macroscopic  physical  parameters.  The  effect  of 
material  gradient  and  FGM  percentage  on  the 
macroscopic  stresses  is  mainly  due  to  the  change  of 
the  distribution  profiles  in  the  middle  parts  of  Figs  8 
and  11,  which  have  lower  magnitude  of  stresses. 

For  the  small  plastic  relaxation  for  local  stress 
concentrations  shown  in  Fig.  13,  this  can  also  be 
related  to  the  fact  that  peak  tensile  stresses  always 
occur  in  metal  grains  where  a  lot  of  ceramic  grains  are 
surrounded.  The  general  stress  state  induced  in  such 
metal  grains  is  mostly  all  around  tension,  i.e.  G\  %  o2 
where  0\  and  c2  are  the  two  in-plane  principal 
stresses.  With  cr,  j= t  o2,  the  maximum  shear  stress 
*^max  ”  (ffi  —  <72)/2  is  therefore  small.  No  matter  how 


low  the  yield  strength  is  in  the  metal  grains,  small 
shear  stresses  can  hardly  produce  any  plastic 
deformation,  and  therefore  are  not  helpful  in  relaxing 
this  kind  of  stress  state. 

4.  DISCUSSION 

A  physically  based  micromechanical  model  is 
developed  to  study  residual  stress  distributions  and 
concentrations  in  the  FGM  sandwiched  between  two 
dissimilar  materials  during  thermal  loadings.  The 
results  obtained  reveal  detailed  information  of 
microstructural  behavior  at  the  grain  size  level,  and 
thus  provide  some  insights  for  optimizing  FGMs  and 
the  control  of  their  failures. 

It  is  stressed  that  the  scale  level  we  are  concerned 
with  in  this  study  is  at  the  grain  size  level  and 


upwards.  Using  the  relatively  coarse  mesh  and 
square  “building  blocks”  is  a  first-order  approxi¬ 
mation  to  the  reality,  and  not  intended  for  subgram 
level  microscopic  features.  For  example,  a  very 
refined  mesh  will  show  very  high  stress  concen¬ 
trations  (in  fact  singularities)  at  the  inter-phase 
corners.  Due  to  the  relative  small  Em.oJE*  rat'° 
(1.77),  the  singularity  is  a  fairly  weak  one,  which 
means  high  stress  regions  will  only  occupy  a  very 
small  volume  fraction.  It  is  shown  that  in  the  case 
of  two  elastic  bonded  quarter  planes,  the  order  of 
the  singularity  is  about  0.05  (1  /»■*.««  0.05)  for  such 
an  elastic  difference,  see  Bogy  [30).  This  kind 
of  singularity  is  thus  expected  to  have  little  (i  e. 
second  order)  effect  on  the  averaged  grain  level 
stress.  However,  if  5-6  vol%  of  the  FGM  grains 
have  some  averaged  tensile  stresses  higher  than 
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Fig.  1 1.  Averaged  in-plane  principal  stress  (j>  ~  (crn  +  <r;;)/ 
2),  distribution  profiles  developed  within  the  FGM  layer 
(40  vol%  FGM,  linear  gradient)  for  both  the  elastic  and  the 
elastoplastic  case. 


Fig.  12.  Macroscopicailv  averaged  in-plane  normal  stress 
developed  within  the  discrete  microstructure  (40  vol% 
FGM,  linear  gradient),  for  both  the  elastic  and  the 
elastoplastic  case. 


700-800  MPa,  only  higher  local  stresses  can  be 
expected  near  those  subgrain  level  microscopic 
features. 

The  grain  level  microscopic  stress  concentrations 
are  found  to  be  quite  high,  of  the  order  of  800  MPa, 
with  only  a  300CC  temperature  drop;  whereas  the 
macroscopic  stresses  are  much  lower  than  the 
microscopic  stresses.  This  suggests  that,  if  high  tensile 
stress  concentration  at  the  grain  size  level  is  the 
failure  initiation  mode  at  this  size  scale,  there  are 
always  micro-fractures  at  the  grain  size  scale  during 
thermal  loading.  Since  the  stress  distribution  profile 
for  the  high  stress  region  is  quite  insensitive  to 
material  gradient  and  FGM  volume  percentage  as 
shown  in  Section  3,  the  above  conclusion  is 
independent  of  those  two  macroscopic  parameters. 
On  the  other  hand,  the  high  stress  region  in  the 
stress  distribution  profile  is  relatively  small  (about 


5-6  vol %  with  stress  p  above  700  MPa).  Whether  the 
small-scale  micro-fractures  at  the  grain  size  level  will 
develop  into  large-scale  fractures  to  cause  the  fatal 
damage  of  the  whole  structure  depends  on  other 
factors  whose  effects  require  further  investigation, 
such  as  stress  redistribution  after  small-scale  failure 
initiation,  loading  history  and  the  grain  boundary 
adhesion  between  the  adjacent  grains.  The  present 
computational  micromechanics  model  can  be  ex¬ 
tended  to  account  for  these  influences. 

For  optimizing  the  microstructure,  since  the  dense 
structure  results  in  high  local  stress  concentrations, 
our  results  suggest  that  for  achieving  higher 
toughness  the  porous  microstructure  should  be 
considered;  and  this  is  consistent  with  the  experimen¬ 
tal  observation  by  Sohda  et  a!.  [11]  on  the  Si-C 
FGM.  Their  experiments  showed  that  the  porous 
microstructure  has  a  much  better  resistance  against 


Fig.  13.  The  6  vol%  averaged  peak  stress  of  p  (APSP)  developed  in  FGM  microstructures  with  different 
material  gradients  and  different  FGM  volume  percentages  for  both  the  elastic  solution  (marked  with  “El" 
on  top)  and  the  elastoplastic  solution  (marked  with  “PI"  on  top). 
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cracking  than  the  dense  microstructure.  Also,  in  a 
recent  review  paper.  Koizumi  and  Niino  [31]  listed 
micropore  as  one  of  the  most  important  material 
constituencies  in  FGM  microstructures.  Regarding 
the  grain  size,  the  extent  of  micro-fracturing  exhibits 
large  sensitivity  to  the  grain  size  in  ceramic 
polvcrystals  subject  to  thermal  loading,  being  more 
severe  in  coarse  grained  ceramics  [32,  33];  therefore, 
a  fine  grain  sized  microstructure  is  suggested  to 
improve  the  FGM's  resistance  against  cracking  and 
delamination. 

Other  microstructural  factors,  such  as  grain  shape, 
grain  size  and  third-phase  particle,  are  beyond  the 
scope  of  the  present  work,  and  are  the  considerations 
of  later  studies.  The  results  shown  here  have  clearly 
demonstrated  the  need  to  consider  microstructural 
details  in  modeling  FGMs  for  determining  their 
mechanical  behavior:  due  to  the  microstructural 
discreteness,  local  residual  stress  concentrations  play 
a  very  important  role  in  failure  initiation. 
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A  Simplified  Method  for 
Calculating  the  Crack-Tip  Field 
of  Functionally  Graded  Materials 
Using  the  Domain  Integral 

A  finite  element  based  method  proposed  for  calculating  stress  intensity  factors  of 
functionally  graded  materials  (\GMs).  We  show  that  the  standard  domain  integral 
is  sufficiently  accurate  when  applied  to  FGMs;  the  nonhomo  gene  ous  term  in  the 
domain  integral  for  nonhomo getfe ous  materials  is  very  small  compared  to  the  first 
term  ( the  standard  domain  integral ).  In  order  to  obtain  it,  the  domain  integral  is 
evaluated  around  the  crack  tip  using  sufficiently  fine  mesh .  We  have  estimated  the 
error  in  neglecting  the  second  term  in  terms  of  the  radius  of  the  domain  for  the 
domain  integration ,  the  material  properties  and  their  gradients .  The  advantage  of 
the  proposed  method  is  that ,  besides  its  accuracy ,  it  does  not  require  the  input  of 
material  gradients ,  derivatives  of  material  properties;  and  existing  finite  element 
codes  can  be  used  for  FGMs  without  much  additional  work.  The  numerical  examples 
show  that  it  is  accurate  and  efficient.  Also ,  a  discussion  on  the  fracture  of  the  FGM 
interlayer  structure  is  given. 


1  Introduction 

The  mechanics  of  functionally  graded  materials  (FGM),  in¬ 
cluding  crack  problems,  have  been  intensively  studied  recently. 
It  has  been  shown  that  for  FGM  crack  problems  the  crack  tip 
has  a  regular  square-root  singularity,  the  stress  and  displacement 
near-tip" fields  are  of  the  same  forms  as  those  for  homogeneous 
materials  (see  Delale  and  Erdogan,  1983,  1988;  Gu  and  Asaro, 
1997a.  b).  So  the  influence  of  material  gradients  at  the  near  tip 
manifests  itself  through  the  stress  intensity  factors.  In  other 
words,  the  stress  intensity  factors  uniquely  characterize  the  near¬ 
tip  field.  Knowing  the  structure  of  the  crack-tip  field,  it  is  im¬ 
portant  to  accurately  calculate  the  stress  intensity  factors  and 
determine  the  effect  of  material  gradients  on  them  for  different 
geometries  and  loadings,  including  those  often-used  specimens. 
Finite  element  analysis  which  can  handle  difficult  material  be¬ 
haviors  and  geometries  as  well  as  various  loadings  provides 
useful  and  the  most  often-used  way  to  solve  mechanical  and 
thermal  problems  including  those  involving  FGMs.  In  this  pa¬ 
per,  we  present  a  simple  and  sufficiently  accurate  finite  element 
method  for  calculating  the  crack-tip  field  for  FGMs,  which  can 
be  easily  incorporated  into  existing  finite  element  codes  and 
commercial  software  packages  without  much  additional  work. 

The  often-used  method  to  calculate  the  crack-tip  field,  stress 
intensity  factors  (elastic  case),  and  energy  release  rate  (elastic 
or  plastic  case),  involves  evaluating  the  7-integral  (Rice,  1968) 
using  the  solved  stress  and  deformation  fields  around  crack  tip. 
For  homogeneous  materials,  this  has  been  an  efficient  way  since 
the  path  independence  of  the  7-integral  allows  us  to  perform 
the  calculation  along  a  path  not  too  close  to  the  tip  so  that  the 
inaccuracy  of  field  variables  at  the  tip  region  due  to  the  singular- 
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ity  can  be  avoided.  Later,  the  domain  integral  method  has  been 
developed  to  perform  the  calculation  of  the  7-integral  (Li,  Shih 
and  Needleman,  1985;  Shih,  Moran  and  Nakamura,  1986; 
Moran  and  Shih,  1987).  The  domain  integral  method  has  been 
shown  to  be  more  efficient  and  more  accurate  than  direct  calcu¬ 
lation  of  the  7-integral,  since  the  domain  integration  comes 
more  naturally  than  the  line  integration  of  the  two-dimensional 
space  and  the  surface  integration  of  the  three  dimensional  space 
in  finite  element  analysis.  Works  along  similar  line  as  the  do¬ 
main  integral  can  be  found  in  the  early  papers  by  Parks  (1974, 
1977),  Hellen  (1975)  and  deLorenzi  (1982),  whose  virtual 
crack  extension  method  is  the  special  case  of  the  domain  inte¬ 
gral.  The  domain  integral  method  has  been  implemented  in 
numerous  programs  to  solve  crack  mechanics,  including  the 
well-known  commercial  package  ABAQUS.  In  this  paper,  we 
use  the  domain  integral  methodology  to  treat  the  FGM  case.  In 
the  non-homogeneous  case,  there  is  an  additional  term  besides 
the  regular  one  due  to  the  variation  of  material  properties.  In 
our  analysis  of  two-dimensional  elastic  crack  problems  of  non- 
homogeneous  materials,  to  capture  the  singularity  and  material 
property  variation,  the  mesh  is  designed  such  that  the  smallest 
elements  at  the  crack  tip  are  very  small,  about  1CT5  times  a 
characteristic  length,  which  is  usually  the  crack  length,  and  even 
much  smaller  for  inelastic  problems.  The  material  variation  is 
achieved  by  using  corresponding  material  properties  at  Gauss 
integration  points  (different  Gauss  points  have  different  proper¬ 
ties).  For  such  mesh  design,  to  perform  the  domain  integration, 
the  domains  can  be  chosen  as  the  circular  regions  formed  by 
the  first  few  rings  of  elements.  In  such  a  situation,  we  show 
that  the  second  term  in  the  domain  integral  for  nonhomogeneity 
is  very  small  compared  to  the  first  term,  the  standard  domain 
integral,  and  may  be  neglected.  Therefore  the  domain  integral 
can  be  calculated  numerically  in  the  same  way  as  that  for  homo¬ 
geneous  materials,  using  the  standard  domain  integral.  From 
the  numerical  point  of  view,  this  allows  us  to  apply  existing 
finite  element  programs  for  homogeneous  materials  to  nonho- 
mogeneous  materials,  avoiding  the  additional  programming 
work. 

The  current  study  is  focused  on  elastic  two-dimensional  and 
three-dimensional  problems.  The  method  may  also  be  extended 
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to  the  nonlinear  material  behavior.  The  numerical  examples 
oiven  include  a  sandwiched  structure  with  a  FGM  interlayer, 
which  illustrates  the  advantage  of  using  FGM  to  reduce  material 
mismatch  between  the  upper  and  lower  layers.  For  such  struc¬ 
tures,  cracks  may  form  at  one  side  and  propagate  to  the  other 
side  through  the  FGM  interlayer  when  the  microscopic  defects 
and  external  loading  are  favorable.  The  crack  in  the  sandwiched 
structure  solved  in’  this  paper  is  along  the  layers’  thickness 
direction  with  the  crack  tip  inside  the  FGM,  and  the  loading 
includes  remote  bending,  three-point  bending,  and  four-point 
bending.  This  kind  of  configuration  may  also  be  good  for  frac¬ 
ture  testing  of  FGMs  since  the  FGMs  are  usually  very'  thin  so 
that  mechanical  testing  may  be  handled  when  bonding  them  to 
two  bulk  materials.  In’  general,  the  solutions  to  the  three-point 
bending  and  four-point  bending  specimens  depend  on  several 
parameters,  including  geometry,  loading,  and  material  variation. 
We  write  them  in  a  compact  form  such  that  the  functionality 
of  each  parameter  may  be  clearly  understood.  These  examples 
show  a  way  to  systemically  present  the  solutions  so  that  they 
can  be  documented  and  are  easy  to  use  in  practice. 

The  paper  consists  of  three  sections.  Besides  this  Introduc¬ 
tion,  Section  2  is  the  discussion  of  the  domain  integral  method 
in  which  we  estimate  the  second  term  due  to  nonhomogeneity 
and  show  that  based  on  the  analysis,  the  term  can  be  neglected. 
Section  3  contains  numerical  examples. 


2  Numerical  Method 

The  crack-tip  stress  field  in  a  FGM  has  a  regular  singularity 
(see  Delale  and  Erdogan,  1983,  1988;  Gu  and  Asaro,  1997a, 
b)  and  the  singular  term  for  plane  problems  is  given  by 
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where  the  angular  functions  and  cr"  are  independent  of  mate¬ 
rial  properties  and  their  variations  and  are  the  same  as  those 
for  homogeneous  materials.  The  displacement  singular  term  is 
siven  by 
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Here,  is  the  shear  modulus  at  the  crack  tip.  The  angular 
functions  uj  and  it"  are  also  independent  of  material  gradients, 
and  are  the  same  as  those  for  homogeneous  materials.  Material 
gradients  only  affect  the  near-tip  fields  through  the  mode  1  and 
mode  II  stress  intensity  factors,  K,  and  K„.  The  energy  release 
rate  is  defined  as 
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which  is  related  to  the  near-tip  field  by 

£0  £o 


(4) 


where  E0  is  Young’s  modulus  at  the  crack  tip.  The  energy 
release  rate  can  be  represented  by  the  following  line  integral: 
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where  IV  is  the  strain  energy  density  and  ns  is  the  outward 
normal  of  the  path  F,  which  starts  from  a  point  on  the  lower 
crack  face  and  ends  at  another  point  on  the  upper  crack  face. 
For  the  homogeneous  case,  the  integral  under  the  limit  is  diver¬ 
gence  free:  therefore,  it  is  path  independent  and  the  limit  is  not 
needed.  In  this  case,  the  7  in  (5)  is  Rice’s  7-integral  (Rice, 
1968).  For  nonhomogeneous  materials,  path  independence  oc¬ 
curs  when  the  crack  is  perpendicular  to  the  material  property 


Fig.  1  A  simply  connected  domain  A,  enclosed  by  the  contour  C  (C  = 
Cl  +  C2  +  C3  +  CA)  near  crack  tip.  The  domain  is  where  the  domain 
integral  is  evaluated. 


variation  direction,  since  in  this  case  the  integral  is  still  diver¬ 
gence  free. 

There  are  usually  several  ways  to  calculate  stress  intensity 
factors  after  the  stress  and  displacement  fields  are  obtained. 
In  the  stress  matching  and  displacement  matching,  the  stress 
intensity  factors  are  obtained  by  extrapolating  from  the  stresses 
or  displacements  ahead  of  the  crack  tip  using  (1)  or  (2).  For 
example,  K,  is  obtained  by  substituting  the  obtained  normal 
stress  ahead  of  the  crack  tip  into  (1).  The  matching  method 
has  the  advantage  that  almost  no  additional  calculation  is  re¬ 
quired  even  in  the  FGM  case,  but  it  requires  a  high  degree  of 
mesh  refinement  and  often  suffers  from  instability  as  the  crack 
tip  is  approached  (see  Anderson,  1995).  Another  way.  the  do¬ 
main  integral  method,  which  is  an  energy  approach  based  on 
the  7-integral  and  which  has  been  proved  to  be  efficient  for 
homogeneous  materials,  is  the  focus  of  our  numerical  study 
here. 

In  the  domain  integral  method,  the  energy  release  rate  7  is 
calculated  through  an  area  integral  in  the  two-dimensional  case 
and  stress  intensity  factors  are  then  obtained  using  (4).  The  area 
integral  approach  provides  much  better  accuracy  than  directly 
evaluating  the  contour  integral  in  (5 ) ,  and  is  easier  to  implement 
numerically.  Early  works  along  the  line  of  the  energy  approach 
were  given  by  Parks  (1974,  1977),  Hellen  (1975),  and  deLore- 
nzi  (1982).  Shih  and  his  co-workers  (see  Li,  Shih,  and  Nee- 
dleman,  1985;  Shih.  Moran,  and  Nakamura,  1986;  Moran  and 
Shih.  1987)  formulated  the  domain  integral  methodology  in  a 
general  way.  For  homogeneous  materials,  it  has  been  applied 
in  above  works  to  elastic  and  plastic  material  responses,  me¬ 
chanical  and  thermal  loadings,  and  two-dimensional  and  three- 
dimensional  spaces.  We  will  discuss  the  application  of  the  do¬ 
main  integral  to  nonhomogeneous  materials.  In  particular,  we 
will  show  that  the  integral  term  representing  the  effect  of  nonho¬ 
mogeneity  may  be  neglected  when  evaluating  the  integration  at 
a  region  close  the  crack  tip;  therefore,  the  standard  domain 
integral  for  homogeneous  materials  gives  sufficient  accuracy. 
We  will  discuss  the  elastic  case;  the  conclusion  may  be  extended 
to  the  power-law  hardening  case,  i.e.,  HRR  singularity  (  Hutch¬ 
inson,  1968;  Rice  and  Rosengren,  1968). 

Consider  an  annular  region  A,  around  the  crack  tip  in  the 
two-dimensional  case,  as  shown  in  Fig.  1.  For  simplicity  in  the 
discussion,  we  consider  that  the  material  variation  is  along  the 
.Y-axis;  and  only  one  of  the  two  material  parameters,  the  Young’s 
modulus,  has  a  gradient  where  the  Poisson’s  ratio  is  taken  as  a 
constant  since  its  variation  is  usually  small  compared  to  the 
former.  The  conclusion  obtained  below  can  be  extended  to  the 
general  material  variation  case.  Both  the  inner  and  outer  bound¬ 
ary  of  the  region  are  sufficiently  close  to  the  crack  tip.  The 
7  given  in  (5)  can  be  written  in  terms  of  the  boundary  integral. 

7  -  (^)  {<JijUiA  -  \Vdij)qnijds  (6) 

where  C  =  C2  +  C3  +  C,  +  CA  is  the  boundary  of  A,:  m.  is  the 
outward  normal  of  Aj\  on  Cj,  nij  —  —Uj,  and  on  C?.  ///.  =  u;. 
and  q  is  a  smooth  function  which  has  the  value  oi  unity  on  (  i 
and  zero  on  CY  Applying  the  divergence  theorem  to  (0)  gives 
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Fig.  2  Finite  element  mesh  of  the  crack-tip  region.  In  our  calculation 
four-node  bilinear  elements  are  used.  The  smallest  element  at  the  tip  is 
10  5  times  a  characteristic  length. 
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Here,  IV  =  VV[E(.v),  c(a\  y)].  The  derivative  of  W  under  the 
second  integral  is  with  respect  to  the  coordinate  .v  in  £(.v). 
Comparing  with  the  homogeneous  case,  the  second  integral  is  an 
additional  term  which  represents  the  effect  of  nonhomogeneity. 

In  numerical  implementation,  the  inner  contour  Cx  is  usually 
taken  as  the  crack  tip,  and  the  outer  boundary  C2  is  taken  to  be 
the  same  as  element  boundaries.  The  function  q  defined  above 
is  an  arbitrary  function  as  long  as  it  gives  the  correct  values  at 
the  boundaries,  Cj  and  C>.  It  was  shown  by  Shih  and  his  co¬ 
workers  in  the  previously  mentioned  papers  that  the  calculated 
J  is  insensitive  to  the  choice  of  q .  The  value  of  it  within  an 
element  may  be  taken  as 


<7  =  X  Ntq, 
1=] 


(8) 


where  Nt-  are  the  shape  functions  of  the  element,  n  is  the  number 
of  nodes  per  element,  and  qf  are  the  nodal  values  of  <?,  which 
are  assigned  in  accordance  with  a  smooth  function  varying  from 
zero  at  the  outer  boundary  to  unity  at  the  crack  tip.  The  deriva¬ 
tive  of  q  with  respect  to  the  coordinate  .v,  is 
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where  T)k  are  the  coordinates  in  the  isoparametric  space.  Evaluat¬ 
ing  the  quantity  under  the  integral  in  ( 7 )  at  the  Gauss  integration 
points,  J  is  obtained  numerically  by 
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Here,  w*p  is  the  weight  function  of  integration,  and  det  (*)  is 
the  determinant  oT  Jacobian  matrix. 

The  mesh  design  for  our  nonhomogeneous  problems  is  a 
standard  mesh  design  for  crack  problems.  The  crack  tip  is  sur¬ 
rounded  by  an  arrangement  of  wedge-shaped  isoparametric  ele¬ 
ments.  The  same  type  of  elements  makes  circular  rings  which 
surround  the  wedge-shaped  elements  at  the  tip  (see  Fig.  2).  In 
this  region,  the  size  of  the  elements  increases  along  the  radial 
direction  according  to  the  exponential  scale  which  gives  the 
unit  aspect  ratio  of  the  elements.  The  smallest  elements  at  the 
crack  tip  are  smaller  or  equal  to  10"5c/,  where  a  is  a  characteris¬ 
tic  length.  Between  the  circular  region  and  the  region  far  away 
from  the  lip  where  the  stresses  vary  regularly,  there  is  a  transi¬ 
tion  zone  in  which  the  element  shape  changes  gradually  from 


the  curved  polygon  to  the  regular  element  shape.  The  geometry 
of  a  typical  mesh  of  this  kind  was  shown  in  the  previously 
mentioned  papers  (also  see  Shih  and  Asaro.  19S9).  Jt  is  noted 
that  near  the  tip  the  mesh  needs  to  be  refined  to, account  for  the 
high-stress  gradients  associated  with  the  singularity;  in  the  FGM 
case,  also  to  account  for  the  material  property  variation. 

Using  this  type  of  mesh,  we  show  here  that  if  we  evaluate 
the  domain  integral  in  the  region  sufficiently  small  around  the 
crack  tip,  the  value  of  the  second  term  in  (7)  involving  the 
derivative  of  VV  is  very  small,  essentially  negligible.  The  domain 
integral  in  practice  can  be  calculated  in  the  region  close  to  the 
crack  tip  (the  circular  domain  consisting  of  the  first  10  or  20 
circular  rings  of  elements  at  the  crack  tip  zone)  as  demonstrated 
in  the  next  section.  In  our  calculations,  as  mentioned  above, 
the  smallest  element  is  in  the  size  of  10~5rt,  where  a  is  the 
characteristic  length.  In  such  a  situation,  the  second  integral  in 
(7)  may  be  estimated  as  follows.  Using  the  above  mesh  design, 
the  first  10  or  20  rings  of  elements  are  arranged  within  the 
radius  10 ~Aa  from  the  crack  tip.  The  weight  functions  for  the 
two  integrals  in  the  expression  (7)  are  q  and  qj.  If  the  pyramid 
shape  for  the  function  q  (Fig.  3(a))  is  used,  its  derivative  with 
respect  to  the  coordinates  is  on  the  order  of  104n"!  considering 
that  the  domain  is  within  a  circle  with  radius  10~4fl.  Then,  the 
weight  functions  of  the  first  and  second  integrals  are  of  the 
orders  104  and  1,  respectively.  Note  that  a']  in  the  derivative 
of  q  has  been  moved  to  the  integrand  of  the  first  integral.  The 
first  integral  is  overweighted  by  its  weight  function  compared 
to  that  of  the  second.  On  the  other  hand,  the  two  integrands  are 
not  likely  to  differ  by  such  a  large  amount  as  that  of  the  weight 
functions,  i.e.,  to  be  of  the  same  order.  This  is  due  to  the  follow¬ 
ing:  (a)  they  both  are  essentially  energy  density  terms  (energy 
density  unit/length)  calculated  using  stress  and  strain  fields; 
(b)  both  are  proportional  to  the  loading,  the  square  of  the  stress 
intensity  factor  K]\  and  (c)  the  first  is  proportional  to  the 
inverse  of  the  modulus  and  the  second  is  proportional  to  the 
derivative  of  the  modulus  divided  by  the  square  of  the  modulus. 
Since  such  a  small  domain  for  the  domain  integral  is  well  within 
the  A-dominance  zone,  the  asymptotic  expressions  (1 )  and  (2) 
are  valid  within  it.  The  E-dominance  zone  for  FGMs  has  been 
examined  in  Gu  and  Asaro  ( 1997b),  where  it  has  been  shown 
that  within  the  distance  of  a  few  percent  of  the  characteristic 
length,  the  difference  of  the  stress  fields  of  the  asymptotic  and 
full  solutions  are  within  a  few  percent.  Substituting  ( 1 )  and  (2) 
into  the  two  integrands,  (b)  and  (c)  can  be  confirmed.  This 
permits  us  to  write 


Fig.  3  The  two  often  used  shapes  of  q  function 
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Fig.  4(b) 

Fig.  4  Two  sandwiched  bending  specimens  with  FGM  interlayer:  (a) 
three-point  bending;  (Jb)  four-point  bending 


[  first  integrand]  = 


[second  integrand]  = 


Kj  MO) 

E0ci  r 

KjEq  MO) 
El  r 


(11) 


In  ( 1 1 ),  Eq  is  the  derivative  of  the  Young  modulus  at  the  crack 
tip:  /i  and/:  are  obtained  from  angular  functions  in  ( 1 )  and  (2), 
and  therefore  they  do  not  have  much  effect  on  the  magnitudes  of 
the  two  integrands.  In  our  discussion,  we  assume  that  there  is 
only  mode  Hoading.  However,  one  can  obtain  the  same  conclu¬ 
sion  for  mixed-mode  presentation  through  similar  steps.  Also, 
Ela  is  roughly  proportional  to  E0  as  can  be  seen  from  the 
following  example.  Consider  the  FGM  interlayer  in  Fig.  4  with 
linear  modulus  variation,  we  have  £0  =  (E2  +  £])/ 2  and  El  a 
=  Eoatyl/h  when  the  crack  tip  is  at  the  middle  of  the  FGM 
interlayer,  where  the  subscripts  denote  the  properties  for  mate¬ 
rial  tt  1  and  #2.  and  C  =  (E2  —  E])I(E2  +  £)).  The  variation  of 
the  multiplier  C.  as  E2/El  varies,  is  small.  So.  the  two  scale 
factors  in  the  above  expression  (11)  would  not  differ  much  as 
long  as  a  and  h  do  not  differ  much,  and  this  is  similar  for  other 
material  variation  forms.  Note  that  if  h  is  the  smallest  length 
compared  to  other  dimensional  lengths,  one  could  chose  the 
characteristic  length  a  to  be  h  so  that  ci!h  =  1.  If  E0  and  E’qci 
are  of  the  same  order,  from  the  above  analysis  we  estimate  that 
the  first  integral  is  10J  times  the  second  integral.  In  general, 
the  difference  between  £0  and  Elo  is  not  significant  at  all. 
compared  to  that  of  the  two  weight  functions  and  therefore  the 
second  term  in  (7)  may  be  neglected.  It  may  also  be  shown 
that  the  conclusion  is  true  by  a  similar  step  if  the  plateau  shape 
for  the  cj  function  (Fig.  3(b))  is  used. 

Suppose  that  rD  is  the  radius  of  domain  where  the  domain 
integral  is  evaluated  and  within  it  the  field  can  be  well  repre¬ 
sented  by  the  singular  field  (1 )  and  (2),  then  from  the  above 
analysis  the  error  of  neelectine  the  second  term  can  be  estimated 


as 


rpEpI  Eo 
1  +  t'pEllEo 


(12) 


The  error  is  very  small  if  we  choose  rD  to  be  sufficiently  small. 
The  simplified  method  has  the  following  advantages:  (a)  it 
gives  the  same  accuracy  for  the  stress  intensity  factors  as  the 
domain  integral  for  homogeneous  materials,  as  we  shall  see  in 


numerical  examples:  (b)  the  input  of  the  derivatives  of  the 
material  properties  (gradients)  is  not  required  such  that  in  the 
numerical  implantation  only  material  properties  need  to  he  as¬ 
signed  to  elements  or  Gauss  integration  points  and  this  is  easy 
to  achieve;  and  (c)  no  additional  calculations  are  required  as 
compared  to  the  homogeneous  case. 

When  the  gradient  of  the  thermal  expansion  coefficient  of 
the  FGM  does  not  vanish  and  thermal  loading  (temperature 
change)  is  applied,  the  additional  term  related  to  the  gradient 
of  the  thermal  expansion  coefficient  is  an  integrand  under  the 
second  integral  in  (7)  which  can  be  written  as 


[thermal  term  under  second  integral]  =  K,T0o: 


,  f:A0) 


\r 


(13) 


where  T0  is  temperature  change  at  the  crack  tip,  a'0  is  the  gradi¬ 
ent  of  thermal  expansion  coefficient  at  the  crack  tip,  and  f\  is 
obtained  from  the  angular  function  in  ( 1 ).  This  term  may  also 
be  neglected  due  to  the  following  reasons.  First,  it  is  under  the 
second  integral  in  (7).  As  analyzed  before,  the  first  integral  is 
overweighted.  Second,  the  Mr  factor  in  first  integrand,  given 
in  (11),  is  much  larger  than  the  \ljr  factor  in  (13),  in  the 
domain  with  radius  10-4rt.  Third,  since  the  stress  solution  only 
depends  on  the  ratio  of  the  muduli  of  the  two  bulk  materials 
for  traction  problems  (as  we  shall  see  in  next  section),  the  two 
moduli  may  be  chosen  in  the  calculation  such  that  1  /£T0  in  ( 1 1 ) 
is  in  a  normal  range.  Having  these,  we  can  estimate  that  the 
magnitude  of  (13)  and  that  of  the  first  expression  of  ( 1 1  )  do 
not  differ  much. 

We  have  gone  a  rigorous  way  to  show  that  the  standard 
domain  integral  can  be  directly  used  for  nonhomogeneous  mate¬ 
rials.  A  simple  way  to  argue  this  is  that  since  the  asymptotic 
expressions  ( 1 )  and  (2)  are  the  same  as  those  for  homogeneous 
materials  with  the  material  properties  being  those  at  the  crack 
tip,  there  exists  a  small  homogeneous  zone  which  may  be  re¬ 
garded  as  the  E-dominan:e  zone  so  that  the  standard  domain 
integral  is  valid  there.  Bui  from  this  simple  way  it  is  impossible 
to  obtain  the  above  error  analysis.  When  the  second  term  is 
neglected,  the  expression  to  numerically  calculate  domain  inte¬ 
gral  becomes 


J 


L  X  j  [(<7yK,.| 

At  p~)  l 


W6ij)q.j]  det 


(  14) 


which  is  the  same  as  that  for  the  homogeneous  case. 

When  the  domain  integral  is  obtained,  the  stress  intensity 
factors  can  be  evaluated  using  (4).  If  there  is  only  mode  I 
stress  intensity  factor  at  the  tip  due  to  the  symmetric  material 
properties,  the  geometry  and  the  loading,  it  can  be  evaluated 
directly  from  (4).  If  it  is  a  mixed-mode  problem,  the  interaction 
energy  release  rate  defined  in  Shih  and  Asaro  (1989)  max  be 
evaluated  instead  of  the  energy  release  rate  in  (4).  Using  the 
interaction  energy  release  rate,  modes  I  and  II  can  be  separated. 


3  Numerical  Results 

We  have  extensively  tested  the  numerical  method  using  many 
crack  geometries  and  loadings.  The  results  have  been  compared 
with  those  obtained  by  other  methods,  such  as  displacement- 
matching  and  singular  integral  equations.  All  showed  the 
method  to  be  accurate,  convergent  to  the  correct  solution.  The 
domain  integral  evaluated  from  the  domains  near  the  tip  is 
stable,  independent  of  the  domain  chosen.  The  following  are 
four  examples  that  illustrate  this. 

The  first  example  is  an  edge-cracked  plate  made  of  a  FGM. 
subject  to  remote  constant  strain.  It  has  been  solved  previously 
by  Erdogan  and  Wu  (1993  i.  using  the  singular  integral  equation 
method.  We  use  it  to  check  the  accuracy  of  our  scheme.  The 
second  and  the  third  are  three-point  bending  and  four-point 
bending  specimens  made  of  sandwiched  structures  with  the  in- 
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terlayers  being  FGMs  (Fig.  4).  The  interlayer  is  a  zone  of 
transition  wherein  the  material  properties  change  smoothly  from 
the  upper  layer  to  the  lower  layer.  The  length  2L  is  assumed 
to  be  sufficiently  large  so  that  it  would  not  affect  the  solutions. 
The  height  of  the  bars  is  2 H,  the  crack  length  a.  and  the  height 
of  the  interlayers  2 h .  The  crack  is  perpendicular  to  the  upper  and 
lower  boundaries,  and  its  tip  is  inside  the  FGM.  The  sandwiched 
structures  can  be  used  to  study  either  the  fracture  of  FGMs  or 
the  interface  behavior  when  the  interlayer  thickness  is  small 
compared  to  those  of  the  two  bulk  layers.  The  geometry  of  the 
first  example  is  the  same  as  Fig.  4,  and  the  only  difference  is 
that  it  is  a  single  piece  of  FGM  for  the  constant  strain  problem. 

The  mesh  design  was  discussed  in  the  previous  section.  The 
four-node  bilinear  elements  are  used  in  the  study.  We  have 
extensively  tested  the  numerical  method  by  changing  the  mate¬ 
rial  properties,  the  loadings,  and  the  specimen  dimensions.  We 
also  have  changed  the  size  of  the  domain  by  changing  the 
number  of  rings  of  elements.  All  of  the  stress  intensity  factors 
evaluated  from  the  domain  integral  have  shown  the  accuracy 
of  the  scheme.  The  convergence  study  results  will  be  provided 
in  tabular  form  later  in  the  section.  The  material  property  varia¬ 
tion  is  achieved  by  using  corresponding  material  properties  at 
Gauss  integration  points  of  each  element.  We  use  the  software 
package  ABAQUS  to  perform  the  calculation,  and  only  the 
user-subroutine  DM  AT  is  required  for  the  material  variation. 
The  ./-integral  is  also  calculated  using  the  standard  domain 
integral  function  provided  in  ABAQUS. 

For  real  FGMs,  the  property  variation  along  the  thickness 
can  be  linear,  exponential,  or  some  other  form.  For  elastic  prob¬ 
lems,  both  Young's  modulus  and  Poisson’s  ratio  vary  with  the 
position  in  general.  It  is  assumed  in  all  examples  in  this  section 
that  the  former  has  the  major  effect  and  the  latter  is  taken  to 
be  constant.  It  is  reasonable  to  do  so  since  the  variation  of  the 
Poisson’s  ratios  is  usually  small  compared  to  that  of  the  moduli. 
For  the  problems  studied  (see  Fig.  4),  the  Young's  modulus  is 
expressed  by  the  following: 

£(y)  =  Ay  +  B 

£(.v)  =  A  exp  (By)  (15) 

where  A  and  B  are  material  constants  which  represent  material 
gradients.  The  origin  of  the  coordinates  is  at  the  center  of  these 
specimens  and  y  is  along  the  thickness  direction.  The  first  ex¬ 
pression  in  (15)  is  a  linear  form,  whereas  the  second  is  an 
exponential  form.  Given  the  moduli  of  the  lower  and  upper 
layer,  E x  and  E:.  the  two  constants  are  expressed  as 


for  the  linear  gradient,  and 

A-m.  «-5;lo.(f)  C7) 

for  the  exponential  variation. 

In  the  first  example,  material  variation  along  the  thickness  is 
taken  to  be  the  exponential  form  in  (15).  The  loading  is  a 
constant  strain  e0  far  away  from  the  crack  at  the  two  ends  which 


gives  rise  to  a  remote  stress  field,  o  —  Oq  exp  (By),  where  cr{] 
=  A(J{  1  -  v2).  The  energy  release  rate  was  calculated  from 
the  domains  formed  by  the  first  20  rings,  using  the  7-integral 
evaluation  function  in  ABAQUS.  The  results  from  the  first  ten 
rings  are  shown  in  Table  I,  given  for  different  Ei/E]  ratios,  and 
the  results  from  the  domains  formed  by  the  remaining  ten  rings 
of  elements  basically  are  the  same  as  those  of  column  9  and  10 
in  the  table.  The  ratio  £;>/£,  is  the  modulus  of  the  upper  bound¬ 
ary  over  that  of  the  lower  boundary.  In  the  calculation  we  choose 
a0  =  1  and  the  crack  length  a  =  1.  From  the  table,  the  conver¬ 
gence  of  the  numerical  method  is  clearly  seen.  When  the  ratio 
is  1,  it  represents  homogeneous  material.  In  the  strong  material 
variation  case,  the  ratio  is  10.  We  see  that  in  both  cases  the 
convergence  behavior  is  the  same.  So,  we  may  conclude  that 
the  convergence  of  using  the  standard  domain  integral  for  FGMs 
is  the  same  as  that  for  homogeneous  materials.  The  domain 
integral  scheme  has  been  proved  to  be  a  useful  one  in  the 
numerical  analysis  of  the  homogeneous  fracture.  Note  that  those 
results  from  the  domains  formed  by  the  first  two  rings  usually 
have  a  relatively  large  error  in  both  the  homogeneous  and  FGM 
case,  due  to  the  inaccuracy  of  the  innermost  elements.  Thus, 
those  results  from  the  domains  formed  by  the  first  two  rings 
may  be  disregarded.  The  stress  intensity  factors  obtained  from 
the  four  cases  in  the  table,  in  which  the  ratio  is  equal  to  0.1, 
0.2,  5,  and  10,  are  4.01, 4.22,  6.49,  and  7.48,  respectively.  These 
results  are  the  same  as  those  obtained  from  the  singular  integral 
equation  method  provided  by  Erdogan  and  Wu  (1993). 

Due  to  the  symmetric  geometries  and  symmetric  material 
properties  with  respect  to  the  crack  line  for  the  two  specimens 
in  Fig.  4,  there  is  only  mode  I  stress  intensity  factor  at  the  crack 
tip  for  the  second  and  third  examples.  Since  the  near-tip  fields 
are  of  the  same  form  as  those  for  homogeneous  materials,  the 
generic  form  of  the  stress  intensity  factor  may  be  written  as 

K,  =  Ta'nY(—,-,-,p)  (>8) 

where  Fisa  representative  stress  magnitude,  a  is  a  characteristic 
length  (can  be  taken  as  the  crack  length)  and  Y  is  a  dimen¬ 
sionless  function  which  is  related  to  the  geometries  of  the  prob¬ 
lems  and  material  properties:  the  ratio  of  the  moduli  and  p.  the 
form  of  material  variation.  There  are  four  independent  variables 
in  the  dimensionless  function  Y.  For  known  material  variation 
p  and  the  thickness  of  the  FGM  h/H,  the  solution  Y  depends 
on  the  modulus  ratio  and  the  position  of  the  crack  tip  a/H ,  and 
may  be  syslemically  presented  by  tables  or  figures.  For  example, 
if  using  tables,  each  table  contains  the  solution  for  given  p  and 
/?///,  where  the  row  represents  E:IE}  and  the  column  represents 
a/H.  If  using  figures,  each  figure  contains  the  solution  for  given 
p  and  /?///,  and  in  the  figure  each  curve  corresponds  to  a  value 
£;>/£,  with  the  jr-axis  being  a/H.  Given  the  representatives  of 
p  and  h/H,  we  construct  the  complete  solution  in  above  ways 
for  others  to  use.  For  homogeneous  materials,  Y  is  only  related 
to  a/H.  It  is  obtained  in  terms  of  a  figure  or  empirical  ex¬ 
pressions  which  are  given  in  the  handbook  by  Tada  et  al. 
(1985). 

Figure  5  shows  the  solution  of  mode  I  stress  intensity  factor 
versus  the  position  of  the  crack  tip  in  the  FGM  for  linear  mate- 


Table  1  Convergence  for  the  Remote  Constant  Strain  Problem 

Energy  release  rate  calculated  from  the  first  ten  rings  [in  units  of  a;n/Ex] 


£,/£, 

1 

-) 

3 

4 

5 

6 

7 

8 

9 

10 

0.1 

42.43 

46.15 

46.21 

46.24 

46.25 

46.26 

46.26 

46.26 

46.28 

46.27 

0.2 

33.26 

36.17 

36.22 

36.24 

36.25 

36.26 

36.26 

36.26 

36.27 

36.27 

] 

20.99 

22.83 

22.86 

22.88 

22.88 

22.89 

22.89 

22.89 

22.90 

22.90 

5 

15.76 

17.14 

17.16 

17.17 

17.18 

17.18 

17.18 

17.18 

17.19 

17.19 

10 

14.78 

16.07 

16.09 

16.10 

16.11 

16.11 

16.11 

16.11 

16.12 

16.12 
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(a-H)/H 

Fig.  5  Stress  intensity  factor  versus  crack-tip  position  in  the  FGM  inter¬ 
layer  for  three-point  bending  with  h/H  -  0.1 


rial  variation  in  the  three-point  bending  specimen,  where  h/H 
=  0.1.  The  geometry  represents  the  case  where  the  interlayer 
of  FGM  is  considerably  thin  compared  to  the  two  bulk  materials. 
It  is  seen  that  the  curves  in  the  figure  are  the  nondimensiona] 
function  Kin  ( 18)  if  the  characteristic  length  is  taken  to  be  H. 
The  solutions  of  this  kind  for  various  h/H  and  p  form  a  com¬ 
plete  solution  for  the  three-point  bending  specimen.  Usually 
tough  materials  such  as  metals  have  a  lower  modulus  than  brittle 
materials  such  as  ceramics.  From  this  figure,  when  the  crack 
tiavels  from  a  tough  side  (the  side  with  smaller  modulus)  to  a 
brittle  side  (the  side  with  larger  modulus)  the  crack-tip  stresses 
increase.  The  energy  release  rate  calculated  from  the  domains 
formed  by  the  first  ten  rings  of  elements  is  listed  in  Table  2  for 
both  linear  and  exponential  material  variations.  The  stable  re¬ 
sults  in  the  table  again  show  the  convergence  of  the  numerical 
scheme.  When  the  toughness  of  the  two  bulk  materials  is  differ¬ 
ent,  it  is  expected  to  vary  along  the  thickness  of  the  FGM  and 
can  be  written  as  T((a  -  H)/H)  in  the  FGM.  Then,  for  stable 
growth  in  the  FGM  interlayer  we  have 


For  unstable  growth,  in  the  second  equation  “<!f  is  replaced 
by  “>/*  Let's  consider  a  special  case  where  the  touahness  is 
constant  across  the  thickness  of  the  FGM.  From  the  figure.  we 
see  in  this  special  case  that  when  material  #2  is  much  softer 
than  material  #1,  E2/E ,  I,  the  crack  growth  is  likelv  to  be 
stable.  This  is  especially  true  when  the  crack  tip  is  close  to 
material  #2.  When  material  #2  is  stiffer  than  material  #].  the 
crack  growth  is  likely  to  be  unstable.  In  general,  if  the  toughness 
varies  with  position  and  the  crack  is  close  to  material  #2  with 


E2/E]  <S  1,  it  is  a  stable  growth  when  the  decrease  of  the  slope 
in  the  figure  overcomes  the  decrease  of  the  toughness.  Figure 
6  represents  the  case  when  h/H  =  0.5  and  other  parameters  are 
the  same  as  Fig.  5.  From  the  twTo  figures  we  see  that  the  trend 
ofthc.se  cur\es  has  a  dramatic  change  as  the  percentage  of  the 
FGM  changes,  and  this  is  especially  true  for  those  curves  with 
E2/E\  <  1.  For  many  crack-tip  positions  in  Fig.  6  the  stress 
intensity  factor  increases  as  the  crack  length  increases.  This 
means  that  if  the  increase  of  the  toughness  at  the  crack  tip  as 
the  crack  length  increases  is  not  as  fast  as  the  stress  intensity 
factor,  it  is  an  unstable  growth  for  the  crack  tip  traveling  in 
those  positions.  In  the  above  discussion  of  the  crack  growth, 
we  have  assumed  that  the  crack  propagates  along  the  original 
direction,  since  these  are  the  cases  where  geometry,  loading, 
and  material  are  symmetric  with  respect  to"  the  crack  line.  A 
first-order  approximation  model,  which  is  based  on  local  homo¬ 
geneity,  has  been  used  to  examine  the  crack  propagating  direc¬ 
tion  for  several  cracked  FGM  geometries  (Gu  and  Asaro. 
1997b).  The  model  also  predicts  that  a  crack  grows  along  its 
original  direction  when  everything  is  symmetric. 

For  the  four-point  bending  specimen  shown  in  Fig.  4.  the 
two  ends  far  away  from  the  crack  line  is  in  a  pure  lending 
state,  where  the  bending  moment  M  =  PI  and  I  is  the  distance 
between  the  applied  force  and  the  support  of  the  beam.  The 
neutral  axis  changes  with  material  properties,  material  variation, 
and  layer  s  thicknesses  and  is  known  when  these  are  given.  Our 
numerical  results  show  very  good  convergence  as  those  listed 
in  Table  1  and  2.  Figures  1(a)  and  (b)  illustrate  the  stress 
intensity  factor  versus  the  interlayer  thickness  for  the  crack  tip 
at  the  center  of  the  beam  and  linear  material  variation,  it  is 
recognized  that  these  curves  are  the  function  Y  in  (18)  if  the 
characteristic  length  is  chosen  as  H .  Compared  to  Fig.  5.  this 
is  another  way  to  present  the  solution:  each  curve  represents  a 
case  for  a  value  E2/Eu  given  p  and  a/H.  In  these  two  figures, 
for  a  fixed  H  and  as  h/H  increases,  the  K}  increases  when  EH 
E\  >  l  and  decreases  when  EHE ,  <  ).  When  h/H  is  zero,  it 
is  the  bimaterial  solution.  The  figure  tells  us  that  although  the 
increase  of  the  percentage  of  the  FGM  in  the  structure  due  to 
the  increase  of  interlayer  thickness  reduces  the  mismatch  be¬ 
tween  the  two  bulk  materials,  it  may  not  reduce  the  stress  inten¬ 
sities  at  the  crack  tip.  The  increase  or  decrease  of  the  stress 
intensities  depends  on  the  crack  position.  As  far  as  the  crack 
propagation  is  concerned,  at  this  point  we  do  not  know  how 
the  toughness  changes  with  the  FGM  percentage  increase  (note 
that  the  increase  of  the  percentage  reduces  the  material  gradi¬ 
ents).  So  we  do  not  know  if  the  increase  of  FGM  percentage 
can  prevent  the  crack  growth.  However,  under  the  special  case 
that  the  toughness  is  constant  for  FGM,  we  may  conclude  that 
the  increase  of  FGM  percentage  may  not  be  good  to  prevent 
crack  growth.  In  another  paper  (Dao  et  a!.,  1997)  we  have 
shown  that  the  increase  of  FGM  percentage  may  not  reduce 
miciostiess  concentration  at  the  grain  size  level  under  thermal 
loading  for  a  perfect  FGM  without  cracking.  The  conclusion  of 
the  microstress  concentration  was  obtained  from  a  statistically 
based  analysis. 

Using  a  similar  analysis,  the  method  can  be  extended  to  the 
three-dimensional  case.  In  three-dimensional  space,  the  domain 
can  be  chosen  such  that  its  boundary  is  a  tube,  which  surrounds 
the  ciack  tip  and  whose  radius  of  the  cross  section  is  sufficient!  v 
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Fig.  6  Stress  intensity  factor  versus  crack-tip  position  in  the  FGM  inter¬ 
layer  for  three-point  bending  with  hlH  =  0.5 


Fig.  8  Penny-shaped  crack  in  a  cylindrical  FGM  solid 
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small  such  that  the  analysis  similar  to  the  two-dimensional  case 
can  be  applied.  Detailed  discussion  of  the  three-dimensional 
standard  domain  integral  can  be  found  in  previously  mentioned 
papers  on  the  domain  integral.  The  reason  for  which  the  simpli¬ 
fied  method  is  true  is  exactly  the  same  as  that  for  the  two- 
dimensional  case,  i.e.,  the  standard  domain  integral  is  over¬ 
weighted  compared  to  the  nonhomogeneous  terms.  Here  we 
will  only  give  an  example  of  a  special  version  in  three-dimen¬ 
sional  space,  an  asymmetric  problem  shown  in  Fig.  8.  The 
detailed  discussion  of  the  general  three-dimensional  problem 
will  be  given  in  a  separate  article  elsewhere.  The  penny-shaped 
crack  (Fig.  8)  is  in  a  cylindrical  solid.  The  radius  of  the  crack 
is  a ,  the  radius  of  the  cross  section  is  R,  and  the  length  of  the 


Fig.  7(6) 

Fig.  7  Stress  intensity  factor  versus  thickness  of  the  FGM  interlayer  for 
four-point  bending 


solid  is  so  large  as  not  to  affect  the  solution.  The  Young's 
modulus  is  considered  to  vary  along  the  radius  direction,  E  ~ 
£(r).  By  the  structure  of  the  crack-tip  field  and  dimensional 
consideration,  the  solution  is  given  in  the  form 

where  E2  is  the  modulus  at  the  outer  boundary  and  £,  is  the 
modulus  at  the  asymmetric  line.  The  solution  to  the  problem 
can  be  easily  documented  compared  to  the  three-layer  bending 
bars,  since  we  have  fewer  dependent  parameters  in  the  nondi- 
mensional  function  Y  here.  For  given  material  variation  p,  the 
solution  can  be  presented  by  a  table  or  a  figure.  Table  3  shows 
the  convergence  for  linear  and  exponential  variations,  where 
E2JE\  =  20  and  ci!R  =  0.5.  The  convergence  is  very  much 
similar  to  the  two-dimensional  case.  Except  for  the  first  one  or 
two  rings,  others  give  the  accurate  solution. 

4  Concluding  Remarks 

We  have  shown  that  the  standard  domain  integral  can  be 
used  to  evaluate  the  crack-tip  field  for  nonhomogeneous  materi¬ 
als.  such  as  FGMs.  The  method  requires  a  sufficiently  fine  mesh 
near  the  crack  tip  as  shown.  However,  the  error  induced  by  the 
method  is  estimated  such  that  one  can  control  the  error  by 
controlling  the  size  of  the  domain  where  the  domain  integral  is 
evaluated.  From  the  numerical  solutions  given  in  the  previous 
section,  we  have  seen  that  the  energy  release  rate  calculated 
from  the  domains  formed  by  the  rings  of  elements  around  the 
crack  tip  in  this  way  is  very  stable  and  accurate.  The  examples 
are  all  in  mode  I  where  both  loading,  geometry,  and  material 
variation  are  symmetric  with  respect  to  the  crack  face,  but  the 
method  can  be  used  to  calculate  the  modes  I  and  II  stress  inten¬ 
sity  factors  for  the  mixed-mode  case  using  a  defined  interaction 
energy  release  rate  by  Shih  and  Asaro  (1989).  The  method  may 
also  be  extended  to  the  nonlinear  case,  such  as  plastic  crack 
problems.  These  all  suggest  that  this  simplified  method,  without 
the  input  of  material  gradients  and  without  many  changes  of 
the  existing  finite  element  code  for  homogeneous  materials, 
may  be  well  suited  for  crack  mechanics  analysis  of  FGMs  where 
the  materials  possess  gradients.  In  the  examples  of  the  sand¬ 
wiched  structure  we  have  presented  the  solution  in  a  compact 
functional  form  which  can  be  used  to  easily  document  a  com¬ 
plete  solution  for  other  study  and  design  purpose.  Finally,  it  is 
noted  that  the  material  variation  is  assumed  to  be  continuous 
across  the  thickness  of  the  FGM  in  this  study.  For  real  FGMs, 
the  material  variation  is  created  by  the  spatial  distribution  of 
one  material  phase  relative  to  the  other.  The  continuous  ap¬ 
proach.  the  proposed  numerical  method  and  the  fracture  behav- 


:e 


Journal  of  Applied  Mechanics 


MARCH  1999,  Vol.  66/107 


Table  3  Convergence  for  Penny-Shaped  Crack  in  a  Cylindrical  Solid*  _ 
Energy  release  rate  calculated  from  the  first  ten  rings  [in  units  of  rWE,] 


*  The  numerical  results  in  the  table  are  for  alR  =  0.5  and  £;/£,  -  20 


ior  analysis  in  the  previous  section,  implies  that  the  particle  size 
of  the  phases  that  make  the  FGMs  is  very  small  compared  to 
the  crack  length  and  other  geometrical  lengths,  and  the  micro- 
structure  of  the  FGMs  is  sufficiently  fine.  In  the  above  condition, 
the  effect  of  the  particle  size  in  the  mechanics  analysis  may  be 
neglected. 
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